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Finite amplitude cellular convection 


By W. V. R. MALKUS and G. VERONIS 
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(Received 22 November 1957) 


SUMMARY 

When a layer of fluid is heated uniformly from below and cooled 
from above, a cellular regime of steady convection is set up at values 
of the Rayleigh number exceeding a critical value. A method is 
presented here to determine the form and amplitude of this 
convection. ‘The non-linear equations describing the fields of 
motion and temperature are expanded in a sequence of inhomo- 
geneous linear equations dependent upon the solutions of the 
linear stability problem. We find that there are an infinite number 
of steady-state finite amplitude solutions (having different 
horizontal plan-forms) which formally satisfy these equations. A 
criterion for ‘relative stability’ is deduced which selects as the 
realized solution that one which has the maximum mean-square 
temperature gradient. Particular conclusions are that for a large 
Prandtl number the amplitude of the convection is determined 
primarily by the distortion of the distribution of mean tempera- 
ture and only secondarily by the self-distortion of the disturbance, 
and that when the Prandtl number is less than unity self-distortion 
plays the dominant role in amplitude determination. ‘The initial heat 
transport due to convection depends linearly on the Rayleigh 
number; the heat transport at higher Rayleigh numbers departs 
only slightly from this linear dependence. Square horizontal 
plan-forms are preferred to hexagonal plan-forms in ordinary 
fluids with symmetric boundary conditions. ‘The proposed 
finite amplitude method is applicable to any model of shear flow 
or convection with a soluble stability problem. 


INTRODUCTION 

This is a study of the non-linear advective processes which determine 
the form and amplitude of cellular convection. We shall pose our problem 
as a formal extension of the conventional linearized stability theory. One 
purpose is to advance a bit closer to the formidable problem of the onset 
of turbulence. 

Linearized stability theory determines those conditions of a known 
steady field of flow which first permit the growth of an infinitesimal 
disturbance. The amplitude of the preferred disturbance or class of 
preferred disturbances is found from the theory to grow exponentially in 
time for values of the external parameters in excess of a critical value. 
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In actuality such disturbances do not grow exponentially without limit, 
but by advecting heat and momentum they alter their own form and the 
distribution of their energy sources to achieve a finite equilibrium amplitude. 
How does this amplitude depend on the external parameters? How much 
heat and momentum are advected? What distortion of the form of the 
disturbance occurs at various amplitudes? What determines the preferred 
state of motion when the stability problem is degenerate, i.e. when more 
than one solution is possible at the point of instability? Under what 
conditions does this finite disturbance itself become unstable? Is the new 
field a similar cellular disturbance or has it the form of an aperiodic 
‘turbulent’ motion ? 

We have chosen to investigate the steady-state, finite-amplitude, vertical 
convection of heat for several reasons. It is perhaps the simplest mani- 
festation of non-linear advection in fluids, both geometrically and in the 
equipment needed for controlled experiments. The stability problem has 
been exactly solved in terms of trigonometric and hyperbolic functions 
(Rayleigh 1916; Pellew & Southwell 1940). Some detailed experimental 
data on the heat transport due to the convection at and beyond the point 
of instability exist (Malkus 1954a) to compare with theoretical deductions. 
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Figure 1. Geometry of the convection problem, 


The physical situation to be studied is shown schematically in figure 1. 
The fluid is contained between two extensive horizontal conducting surfaces 
distance d apart. ‘The upper surface is held at the temperature 7,, and the 
lower surface at the higher temperature 7',. In the absence of motion the 
temperature distribution in the fluid is determined solely by the thermal 
conductivity and ‘s indicated by the heavy line connecting 7,, and T,. 
In 1916 Lord Rayleigh deduced the criterion for marginal stability in such 
a fluid layer. ‘This criterion involves what has come to be called the 
Rayleigh number 


wo as 
A =— Pas d 
KV 


where « is the linear coefficient of expansion of the fluid, g is the magnitude 
of the gravitational acceleration, 8,,=(7,,—7,)/d is the mean negative 
temperature gradient, « is the thermometric conductivity and v the 
kinematic viscosity of the fluid. When A exceeds a critical value convection 
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occurs as a regular cellular pattern. The theory does not distinguish 
between rectangular, triangular or hexagonal horizontal plan-forms for these 
cells. All could occur at the same value of A with identical exponential 
growth rates. Another undetermined quantity in the linear theory is the 
sign of vertical motion in the centre of the hexagons. 

The observations recorded in the earlier experimental paper (Malkus 
19544) establish that for values of \ up to ten times the critical value A, a 
steady cellular convection exists in the fluid. The heat transport due to 
this convection shows a remarkably linear increase from zero at the critical A 
to a value greater than that due to the conduction alone at A=10A,. This 
is the range of A to be studied in this paper. 

At A= 10A, a new instability occurs in the fluid producing disordered 
aperiodic motions, quasi-cellular in appearance. ‘This has been interpreted 
as the onset of some type of turbulence. However, the heat transport 
once again increases linearly with A, but with a steeper slope than in the 
cellular region. Further discrete changes which steepen the slope of the 
heat transport curve appear at higher values of A and appear to be associated 
with further instabilities. Each transition leads to a more intense and 
apparently more disordered field of motion. A preliminary statistical 
study of the convection at these high values of A has been made in an earlier 
paper (Malkus 1954b). Several results of the present study of cellular 
convection which relate to the aperiodic convection at large A will be discussed 
in the conclusion. 

The dashed curve connecting Ty and 77, in figure 1 is an example of the 
distribution of horizontally-averaged temperature when cellular convection 
occurs. ‘The boundary conditions prevent convective heat transport at 
the boundaries. Hence the additional heat transport is permitted by, 
increased mean temperature gradients at the boundary and decreased 
gradients in the mid-regions of the fluid. One task we now undertake is 
to determine the magnitude and shape of this distortion of the mean 
temperature field. 

In the first section we describe the iterative method of solution for the 
finite amplitude fields of temperature and velocity. This is followed by 
a detailed treatment of the simplest case, that of two-dimensional convection 
(roll cell) for free surface boundary conditions. We were able to carry 
this analysis to ‘sixth order’ and investigate many of the ‘distortions’ 
of the initial convective disturbance. In§3 we study the first finite amplitude 
effects of cellular convection with rectangular and hexagonal plan-forms 
(three-dimensional cells). Up to this point a multiply-infinite set of 
steady-state finite amplitude solutions has been generated, and we turn 
to the problem of their stability. The ‘relative stability’ argument of §4 
resolves the indeterminancy of the steady-state equations by providing a 
criterion which selects the realizable solutions. In $5 we discuss the finite 
amplitude solutions when the fluid has rigid boundaries and in conclusion 
describe certain approximation techniques to extend the useful range of 
the formal analysis. 

P2 
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1. BAsiC EQUATIONS AND THE METHOD OF SOLUTION 
‘The equations which have been used by past authors to investigate 
convective instability are attributed to Boussinesq (1903, p. 173). They 
are based on the assumptions that the viscosity and thermal conductivity 
of the fluid can be treated as constants and that variations in the initial 
density field are important only in the buoyancy-force term in the equations 
of motion. Jeffreys (1930) has shown that the resulting disturbance 
equation is valid in compressible fluids if the temperature gradient minus 
the adiabatic lapse rate is used rather than the actual temperature gradient 
and if the total variations in density in the fluid are very small compared 
to the mean density. In liquids these equations for the local conservation 
of heat, momentum, and mass respectively are 


(a/at—nV2)T=—-V.VT, a) 
(e/et —vV2)V = —V.VV— VP/p,) + y(T — Tok, (1.2) 
v.V—0, (1.3) 


where V is the vector velocity of the fluid, T is the actual temperature, 
T,, is a standard temperature and p, a standard density. P = P—gz where 
P is the pressure, y = ag, and k is a unit vector in the z-direction. The 
non-linear advective terms V.VT and V.VV are to be retained in this analysis 
and are the principal objects of interest. However, in applying these 
equations to finite amplitude processes it must be borne in mind that the 
maximum density fluctuation from the mean must be small. Several 
interesting geophysical and astrophysical applications of this work require 
reappraisal of the assumptions and the inclusion of additional terms in the 
equations. 

To clarify the role of the non-linear terms they will be divided into terms 
which are finite when averaged over a horizontal plane and which therefore 
modify the average field, and into terms of zero average. 

Denoting an average over a horizontal plane by a bar above the quantity, 
the heat flux equation (1.1) becomes 


~ + = (xB) =— = (WT) (1.4) 
where 8 “T cz, T=T—T, and W is the vertical component of V. 
Using (1.4), (1.1) may be rewritten as 
(c/ct—KV?)T=BW-A, (1.5) 
where h=V.VT-o0(WT)/ez, h = 0. 
From (1.2) and (1.5) two useful relations can be generated which describe 
the gross energetics of any flow. Multiplying (1.5) by y7’8,,, (1.2) by V, 
ind averaging the resulting equations over the entire fluid, one obtains 


ly cm C sae.) YK (SP eT 6) 

5 Tm = We A r) - & (VENT (1.6) 

and LC ee TR ? 7 7 
(V.V),, = 7(WT),,—v(VV;.VV)) (1.7) 
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where the subscript m denotes the mean value over a vertical line joining 
the two horizontal boundaries, and V; denotes the 7-component of the 
velocity. The averages over the triple-product terms generated by the 
zero average non-linear terms and the average over the work term V.VP 
must all vanish since they represent conservative advections within the 
fluid system (implying that (AV.VA),,=0 when V.V=0 and 4 is any 
fluctuating scalar field). 

Equation (1.7) is the ‘ power integral’ of the motion. It states that the 
time rate of change of kinetic energy per unit mass is equal to the rate of 
release of potential energy by the convection minus the rate of dissipation 
of kinetic energy by the viscous stresses. Equation (1.6) has been written 
as a ‘power integral’ to parallel (1.7). However (1.6) is more correctly 
interpreted as either the balance equation for the mean square of fluctuations 
of the internal energy, or, equivalently, as an entropy balance equation. 
The particular value of these non-linear integral statements is in the 
determination of amplitude when the form of the motion is known or 
adequately approximated (e.g. Meksyn & Stuart 1951). We shall consider 
them later in this paper. 

Since we will pivot our analysis about the linearized stability problem, 
it is convenient to eliminate all but one of the dependent variables appearing 
as linear terms in (1.2), (1.3) and (1.5). Therefore, we cross differentiate 
to remove the pressure and the linear terms involving the two horizontal 
velocity components U and V in (1.2) and (1.3). The resulting relation 
between IV’ and T (plus non-linear terms) is 


(c/ct—vV?)V?W = yViT+L(M), (1.8) 
a2A 92] 
where L(M) = aM. + os, —V? M., 
Oxdz yes 
M,=V.VU, M,=V.VV, M, = V.VW, 
and Vi = 0?/dx? + 0?/dy" 


The two equations which relate the linear components U and V to W are 




















(g--P)(U+ )--Setee 
ct OxX0zs oy* cxcy 

é , Ww o*7M 02M 

=. «eee Urs = =. e. 1.10 
( a ) ( : oy ==) Oxody Ox? ee 


The linear parts of (1.8), (1.9) and (1.10) are just those found in the analysis 
of convective instability. Cross differentiating (1.8) and (1.5) produces 
the usual sixth-order equation in the one ‘linear’ variable W: 
< —xV?)( = —vV?)V?W—yBV?2 W = —yV2h+(= —«V?2)L(M). (1.11) 
ct ct ot ; 
All terms on the right of (1.11) are non-linear terms with zero average. 
On the left are linear terms except for the important horizontal-average 
temperature gradient which depends upon WT. To establish this latter 
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relation we first note that the time-independent integral of (1.4) for fixed 
boundary temperatures 1s 

«B+WT = H = xB, +(WT)ms (1.12) 
where H is the constant heat flux between the horizontal surfaces due to 
both conduction («8) and convection (WT). Then from the general 
integral of (1.4) one may write 














B= 1+ 4 (WP WT +o(—) ‘ (1.13) 


where G is that part of 8 which vanishes when the time variations of WT 
vanish. G is perhaps important in the study of the onset of the aperiodic 
‘turbulent’ convection, but since it will soon disappear from this problem 
its explicit form need not be given here. 

At this point it is of value to non-dimensionalize the physical quantities 
appearing in the equations. ‘This will simplify the mathematical manipula- 
tions and produce the pertinent physical a, in the problem. Let 
V=(x/d)V’, T= (ve/yd*)T’, (x, 9,3)= d(x’,y’,2’), t= (d?/«)t’. (1.14) 
Having made this change all primes will now be dropped. Jn the remainder 
of this work all unprimed quantities are non-dimensional unless it ts otherwise 


stated. Using (1.13) and (1.14), equations (1.6), (1.8) and (1.11) become 
. anes G)n}—(VT.VT) 
(o— d/ot — V2) V?2W = V? 7101 (M) (1.16) 

(0/at — V2)(o-! d/at — V2)V2W i ns WT+GiR2W 
V?h+oa-1(e/et—V*?)L(M)_ (1.17) 


1a(T?),, at = (WT), +{(WT, — (WT? 


mi 


where o = v/« is the Prandtl number, and J is the previously defined Rayleigh 
number. ‘The non-dimensional forms of equations (1.9) and (1.10) will 
be used later in the analysis but are not needed immediately. 

In the usual studies of disturbances of the state of steady conduction, 
one seeks solutions of the form 


W=eW,+2W,+8We+..., 


{ (1.18) 
T ~€8,40T. +74. 


with similar expansions for the other variables, where € is a constant 
parameter. It is then required that expressions (1.18) satisfy the complete 
equations of motion for all values of € less than some maximum e. The 
coefficients of each power of « generated by substituting the expressions 
(1.18) into the equations of motion must vanish individually and the resulting 
series for each of the variables must converge if relations (1.18) are to 
represent a satisfactory solution to the problem. Stability theory is 
concerned with the solutions of the first-order equation only. For first-order 
solutions to be complete one must require that « be proportional to the 
amplitude of the disturbance and that this amplitude be infinitesimal. In 
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this finite amplitude study we will find that a similar identification of « and 
amplitude is necessary if (1.18) is to represent a solution*. 

An important point in this expansion remains to be discussed. When 
the actual forms for the W; and 7; are substituted into (1.15), one finds an 
explicit relation between e and A. Thus 

A= Aypted, +e7Ag+..., (1.19) 
where the numbers ), are integrals of functions of W,and T;.. Sucha relation 
between the amplitude and the physical parameter is to be expected. One 
might have preferred the amplitude to be given as a power series in A since A 
is actually the variable which can be controlled. However, (1.19) is a natural 
consequence of the expansion (1.18). Explicit determination of the A; is 
discussed in the following paragraphs. 

The sequence of linear inhomogeneous equations in powers of € generated 
by inserting equations (1.18) and (1.19) into equation (1.17) is 
L(W,) = (a/et— V2)(o- a /at — V2)V21V, —AV7 Wy = 0, 

L(W,) =r, V2 Wo — V2 Nog +o71(0/0t — V?) Loo; 


~W.T,+GyjV2W,- | 
| 
| 





f(W,) = Ay VW, t X53 ViW,+ (WW To) m 
j Vi(Ao u hy) ne 1(0/dt — V?)(Lo1 ¥ Lio); p (1.20) 


nxt 2 
54,,V2W,+ 5 
0 l 





II 


¢(W;) i-n V CW, Tn — W,, T+ Gh Vi Win sty t 


ni—n+lL) 


tS {-V2hy inant o71[0/0t — V2IL 


where the subscripts n, / on the non-linear terms h, L and G mean that 
W,, or V, and T, are to be substituted for W or V and T in these terms. 
The 7; can be determined from the IV; by the auxiliary equations derived 
from (1.16) 

V2 T, = (o—1 d/at — V2)W, 


oO 

it 5 (1.21) 
Vi T; = (0-1 d/ot — V2) V?W,-o7 ¥ Las naps | 
0 
with similar expansions for U,; and V;. Each of the W; must satisfy the 
boundary conditions for the vertical velocity, each of the 7; the boundary 
conditions for the temperature fluctuations. The various boundary 
conditions of the problem will be discussed in the next section. 

The first of equations (1.20) is the classical Rayleigh stability equation, and 
¥ can be termed the linear constant-coefficient Rayleigh operator. Assume 
tentatively that it is solved for W, (normalized) and A,, and that from the 
auxiliary equations 7, Uy, Vy are also known. The second of equations 
(1.20) is a linear inhomogeneous equation for W,. The inhomogeneous 
terms depend on the known forms of W,, 7), U,, Vj and the unknown 

*V.S. Sorokin (1954) attempted to solve the second order equations by assuming 


that € ~(\/A—v/Ap)'*. His neglect of the important terms containing //f,,, in the 
basic equations prevented him from obtaining any quantitative finite-amplitude results. 


n 
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number A,. ‘The solution of this equation for W, will be the sum of a 
homogeneous solution plus the particular integral generated by these 
inhomogeneous terms. ‘Three difficulties arise at this point. (1) The 
part of the homogeneous term which has the form of W, will produce a 
secular (non-periodic) response in the particular integral—a solution which 
has no physical counterpart; (2) that part of the homogeneous term which 
satisfies the boundary conditions will contain an arbitrary constant; (3) a 
part of the inhomogeneous term may satisfy “(W) = 0 but not the boundary 
conditions (i.e. it will differ from W,) and will also give rise to a secular 
response. ‘The last occurs nowhere in our analysis because the inhomo- 
geneous terms themselves always satisfy the boundary condition. ‘The first 
two complications require that a method be adopted which resolves the 
indeterminacy and which makes the e-equations a soluble iterative sequence. 

We have found only one such method*. It is to require (1) that the d, 
be evaluated so as to eliminate the ‘ resonant’ inhomogeneous terms, for 
example 





= (W, V2 hoo) mn — 7 { Wo(0/et —V?)Lootm> (1.22) 
and (2) that (W,W),, = , i.e. that « be proportional to the amplitude of that 
part of }V’ which has the form Wy. ‘Then from (1.18) all W; (¢ > 0) must 
be orthogonal to MW, though not necessarily orthogonal to each other. The 
latter requirement removes the arbitrary constants of the homogeneous 
parts of the W, equations. 

To eliminate the resonant terms in the second-order equation, we must 
have 


(W,,) 


—1, (Wy Vi Wy) - 
— } Wl Wy To) m ae W,, Ty Tr Goo] Vi Wom ” 
[Wo Vi(hor + 410) Im — 9 [ Wo(0/ et — V*)(L 9, + Ly0) | m- 
(1.23) 


Higher A; generated from (1.20) by the first requirement above are 


mn 














u“ 


oset es eR 
As = (Wo St Won § = DAR ( Wo Vi Win + (Wo Vi bn in tm — 


i U0 





-o"[(W,(0/et—V2)L,, «pln - 








> K(W,, Tim — Wr T+ Gat Wo Vi Weensisolme (1-24) 
7 } 


A formal technique for determining the amplitude and distortion of 
a disturbance as a function of the external parameters is now complete. 
Can one determine the range of A in which a series solution terminated 
after a term e€' is a good approximation to the whole series for W? The 
best we can do to test the validity of a limited solution is to compare the 


* Recently we have learned that an analogous method was proposed by Lindstedt 
(1883) to obviate similar difficulties which arise in problems of celestial mechanics. 
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amplitude of W; to «, the amplitude of W). If this ratio exceeds 0-1, for 
example, we can be fairly sure that an additional function W;,, must be 
computed to keep errors in distortion below 10°,. However, the important 
energy and heat transport terms are proportional to products of W and T. 
Hence the error in these terms will be much smaller than the error of Wand 7 
themselves. Examples of determination of the errors in limited solutions 
are given in the following section. 

The proposed technique cannot answer all the explicit questions raised 
inthe Introduction. Ifasolution tothe stability problem, Wo, isindependent 
of time, then none of the equations (1.20) can lead to time dependence. 
If there are many solutions to the stability problem, then the equations (1.20) 
give as many answers. ‘l’o determine which of these answers is realized in 
an experiment we must separately investigate their ‘ relative stability ’. 

Therefore, we do not present this technique as a unique statement of 
the finite amplitude problem. Several alternative approaches which also 
answer some but not all of our questions are discussed in the Conclusion. 

The following section discusses the consequences of this e-expansion 
with generating functions based on the simplest solution of the Rayleigh 


problem. 


’ 


2. ANALYSIS FOR THE TWO-DIMENSIONAL CASE WITH FREE SURFACES 
The first step in a determination of the effects of finite amplitude in 
equation (1.20) is the resolution of the stability problem. Pellew & 
Southwell (1940) have made a most comprehensive study, which will 
be only briefly outlined here. They investigate the first of equations (1.20), 
#(W,) = 0, to determine W, and A, for three different sets of boundary 
conditions on the velocities. In every case the boundaries are perfect 
heat conductors, that is, the fluctuation temperature vanishes there. The 
conditions at a free boundary are 
W = PWios* = T = 0, (2.1) 
and are a consequence of the divergence relation (1.3) when the boundary 
cannot support a stress. The conditions at a rigid boundary are 
W = dW/dz = T = 0, (2.2) 
and are a consequence of the divergence relation when all velocity components 
must vanish at the boundary. 
The assumption made in solving the Rayleigh problem is that the field 
of motion is separable, 1.e. 
W, = $(x, y)F(z)G(t), ) 
and Vi d(x, y) = —a®nx* d(x,y), J 


4 


(2.3) 


where the separation parameter « is the effective wave-number of the 
disturbance in the horizontal plane. (In what follows we treat (x,y) as a 
normalized function.) Pellew & Southwell have shown that this assumption 
is justified only if the horizontal plan-form of the motion consists of regular 
‘close-packed cells’ at whose lateral boundaries the normal derivatives of 
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the temperature and velocity fluctuations vanish. ‘The plan-forms which 
satisfy these requirements are the two-dimensional ‘roll’, the general 
rectangle, the equilateral triangle and the hexagon. Hence a characteristic 
value A, will have associated with it an infinite number of possible 
characteristic functions, all with the same «. Pellew & Southwell correctly 
suggest that this degeneracy in the first-order theory is removed by ‘ higher- 
order’ effects, though they were apparently under the impression that the 
hexagonal plan-form was the preferred motion in all experiments. 
Another important contribution of Pellew & Southwell to this first-order 
problem is their proof that all oscillatory solutions decay. Hence maintained 
convection is initiated as a steady motion and “‘limiting conditions of 
stability are in fact obtained when all time variations (in 4(W,) = 0) 
are made zero’’ 
Therefore, with the use of (2.3), the equation #(W,) = 0 becomes 
(62/02? — a2?) W, +A, a?n?W, = 0, (2.4) 
which is a linear sixth-order equation with constant coefficients. Its 
general solution 1s 
W, = 4(x,v) & (A;cosh 2y,3 + B;sinh 2p; 3), (2.5) 
7=1 
where A; and B; are arbitrary constants and 
4? = 2? 2?{1 — (Ay/x?2?)!3w;}, (2.6) 
where the w, are the three cube-roots of unity. 
Application of the conditions for two free boundaries to (2.5) and (2.6) 
generates the characteristic functions 


W, = A (x, y)sin nzz, (2.7) 
where m is an arbitrary integer, and the characteristic values are 
Ny, = 74 (n? + x?)8/a2, (2.8) 


The lowest value of Ay,, is 2774/4, and it occurs for x7 = } andn=1. The 

case n = | is our primary concern, since, before the modes corresponding 

to higher values of m can occur, the mode corresponding to n = 1 has 

grown to finite amplitude, markedly altering the basic temperature field. 
For two rigid boundaries the characteristic function becomes 


W, = 6(x,y) & A;cosh 2u,(z— 3), (2.9) 
l 
where 23 = +1 at the boundaries and where A, and A, are related to 4, 
by the conditions 
> A;cosh p; = 0, > »; A;sinh p; = 0. (2.10) 


1 1 
The yu; are found from (2.6) once A, and x are known. The relation between 
Ay and « has been determined numerically by Pellew & Southwell (1940, 
p. 377). ‘The lowest value of Ay is 1707-8 for 2? = 3-13/7. 


For one rigid and one free boundary the first characteristic function 
is the second, or ‘odd’ mode of (2.9) with 2z = 1 at the rigid boundary 
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and 2: = 0 at the free boundary. Pellew & Southwell found the minimum 
eigenvalue to be A, = 1100-65 for «2 = 2-68/z. 

The primary gonclusion drawn from the first-order equations is that 
for A < A, all infinitesimal disturbances to the purely conductive state decay 
in time. What will be the equilibrium amplitude of the preferred 
disturbances for A > Ay? 


Finite amplitude solution for rolls 

In the remainder of this section we will study the solutions generated 
by (1.20) for the very simplest case, that of the two-dimensional roll with 
two free boundaries. Many, but not all, of the properties of the other 
solutions to this problem are exhibited in this case. However, it must be 
borne in mind that neither the two free boundaries nor the roll can be 
realized in practice. 

A complete solution to the first-order equations for strip plan-form 
and two free boundaries is 


W, = 2cos7axsin 73, 
Ty) = (1+ «?)?(27?/a?) cos 7ax sin 72, 
U, = —(2/«)sin 7x cos zz, Vi = 0, | 
Ny = 4(1 + a?)3/x?, 
from (1.20), (1.21), (1.9), (1:10) and (2.1). 
As noted previously the lowest value of Ay is 2771/4 and occurs when 
z* = 3, Before solving the second-order equation for W,, the quantities 
Nog, Log and A, must be determined from (2.11). Now 


hey V¥_. V7,— — (W, le) = | pes 2az— < (1sin? zs | = 0, 


J ~ (V,. VU,) + ~~ (Vp. V¥e)— V2V_. Vo = 0. 
zs CYOR 


I 





i= = 
ONC 


Hence from (1.22) A, = 0. We shall find that A, = 0 for all symmetric 
solutions to the first-order equations. However, the conclusion that 
hoy = Ly) = 9 is true only for rolls and only with two free boundaries. 
This fortuitous vanishing of all zero-average advection by tbe first-order 
solutions leads to “/(W,) = 0; hence W, = U, = 7, = 0 from (1.20) and 
(1.21). Therefore hy, = Ay = Lo, = Ly) = 0, and the third-order equation 
of (1.20) becomes 
L(W,.) = rA. V2? Wo +i (Wo To)m — Wo To Vi Wo, (2.12) 

where A, is given by 

No( Wo Vi Wo) mn = —[6(Wo Tom — Wo To} Wo Vi Wolan (2.13) 
from (1.23). From (2.11) and (2.13) 

















236 W. V. R. Malkus and G. Veronts 


Therefore (2.12) becomes 
V8W,—A, Vi W, = —7*(1+2?)? cos 7ax sin 372. (2.15) 


Before investigating the distortions of the roll form generated by (2.15), 
the implications of our first finite-amplitude results, equation (2.14), will 
be studied. From (1.19), (1.18) and (2.11), to this third-order approxima- 
tion, e? = (A—A,)/A, and 


(WT), = &(W, T,),, = 2(A—Ap).- (2.16) 
Also, from (1.13) in its non-dimensional form, 
B/B,, = 1+2(1—A,/A)cos 27. (2.17) 


In (2.17) the distortion of the mean temperature gradient £/8,,, is determined 
by products of zero-order functions only and requires that for A > 2A, 
the gradient should become negative in the mid-regions of the fluid. 
Equation (2.16) asserts a linear relation between the convective heat 
transport and A. (Note that this relation is independent of o.) ‘That the 
observed heat transport shows just such a linear relation from Ay to the 
second transition at 10Ay suggests that this third approximation provides 
an adequate description of the field of motion throughout the entire 
range of steady convection. However, this is not the case as we shall now 
establish by the investigation of higher order terms. 


Ay-approximation 
The solution of (2.15) for W4 is 


W, = Cy cos wax sin 3x2, (2.18) 
where 


Cy = 
Hence from (1.21) and (1.9) 





T, = C,cosmaxsin3nz, U, = —(3/x)Cy sin waxcos3mz, (2.19) 
where 
Cy = (x/afP(S* +a? fC. 

Therefore the first distortion of a finite-amplitude roll does not affect the 
horizontal plan-form, but only the vertical structure. ‘The addition of the 
term sin 37z to the first-order solution increases the amplitude of W and T 
near the boundaries in keeping with the increased gradients near the 
boundaries (see equation (2.17)). 

This rather straightforward iterative analysis will now be extended to 
determine As, Ay, A; and Ag and the accompanying W,, T;,, U;. 

From (2.11) and (2.19) one finds that 


Vi (hoz + Aap) = Cy(x)cos 27%x(2 sin 273 + sin 472) 


and = a 7 V2(Loo + Lop) = 0 !C,(«)cos 27ax(2 sin 272+3sin47z). (2.20) 
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Since W, is orthogonal to the expressions (2.20), and since W, 7, = W,T, =90, 
then A, = 0. Therefore 
¥(W3) = cos 27ax[C,(x, o)sin 273 + C,(x, o)sin 473]. (2.21) 
IV, and 7, are of the same form as #(W,), but with different coefficients 
C(x, 0), and 
U, = sin 27ax[C;(«, «)cos 27z + C,(a, o)cos 472]. 
In order to determine A, one must construct Jigs, Azo, Lo; and Ly from (2.11) 
and (2.21). ‘This leads to 
Vi (ho; + hgo) = cos 7ax(D, sin 372 + Dy sin 522) + 
+ cos 37ax(D,sinzz+D,sin37z+ D;sin5zz), (2.22) 
with an identical form for V?(Zo3;+ 9) differing only in the coefficients 
D(z,c). Therefore 


[W, Vi (Ags + hzq)|n = [Wo V?(Lo3 + £30) ]m = 9 (2.23) 








and 





[(W, T.+ WT ))Wo Vi Wolmn + (Wo Ty Wo Vi W2) 
(17, V3 Wy), 
= —}Cp+2nX(1 +222 2Cy. (2.24) 


m 





Since A, is an integral of zero-order and second-order functions only, it is 
not dependent upon o. Hence, to the fifth order, o does not influence the 
amplitude of the convection. Indeed the largest part of the o-dependent 
inhomogeneous term in the equation for W,, namely, o-!V?(Lo3 + Lao), 
is smaller than the corresponding term in Vi(hp3+/y9) by a factor 
x? [2n?(1+)’o], = 10-*/o. Hence in the determination of W, and A, 
these higher-order momentum advection terms will be neglected. 


\y-ap proximation 
Proceeding as above one finds that 
W, = cos 7x(F, sin 372+ FE, sin 572) + 
+ cos 37x(E, sinzz+ E,sin3az+E;sin57z), (2.25) 


with a similar form for 7, where E = E(z,c). 





As for previous 4; with odd 7, A, = 0 since W, W, = W, T; = W3T, = 0 
and [W3(Aoy +Ao2)],, = 9. Therefore 
As = (WoViW),, (WoT + WiTo)(WoVi Wo) + (Wey + WoTs) WoVi We + 
+ WyT) WoV? W,— (a/7)2(1 + @?)-2[ 7, V2 (hog + Ago) + 
+ T,V=(Ioo + hoo) ]} (2.26) 


Sm? 








where use has been made of the facts that (Th),,=0 and 7, = N,W, 
in order to simplify the last term. 

Here for the first time the Prandtl number can have an effect upon 
amplitude, since the coefficients in W; and W, are functions of o. 
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We are now able to investigate the ranges of A in which these various 
solutions W,, W;, W, first become significant contributors to the field of 
motion. 

For the initially unstable disturbance A, = 277'/4 and «#2 = 3 and 
from (2.14), (2.18), (2.19) and (2.24) A, = 97?/4 and A, = —0-1248 = —-}. 
Then, to this A,-approximation, from (1.19) 


9 1 § Ae 6 . (A — Ao) oi 7 | 
Rt at i PAO EOE be " yp ae 
i ae | %) i J oa 


Equation (2.27) tells one that e? becomes imaginary when 
A “_ Ny = = AS Ay = 13A,. (2.28) 


Hence the range of A in which this A,-approximation is valid is certainly 
less than 1:5A), while the range in — the previous \,-approximation is 
valid must be smaller still, When A—Ay = 1-5Ag, e? = — SA,/A,y = 2(A—Ap)/Ad 
from (2.27) and (2.28), and e? is twice as great as the amplitude predicted in 
(2.16) for the A,-approximation. ‘Therefore the hope that the A,-approxima- 
tion would be valid throughout the laminar convection range is unjustified. 
It is now clear that the distortion of the initial disturbance must play a 
significant part in determining the observed heat transport. 

The A,-computation from (2.26) though tedious, has proved of 
considerable value in explaining both the observed relation between 
A and (WT),, and the nature of the e-expansion. Since A, is a function 
of o, the various coefficients, C(x,0), D(«,o), E(x,o), entering into it 
have been computed for «? = } at three values of o. For o = 0-8, a value 
appropriate for a gas, Ag = 2:°31x10-%. For o = 8, a value appropriate 
for water, A, = 2°30 x 10-3. For o= o, A, = 2:°29x 10-3. It would be 
necessary to include the small neglected o-term in W, in order to decide 
whether the effect of o on heat transport is just very small or whether it 
vanishes completely. We will return to the effects of o in the following 
section. Here we will study amplitudes and distortions when A, = 2:30 x 10-% 
corresponding to the value of o for water. 

In order to determine ¢€ as a function of A, to the A,-approximation, 
from (1.19) one must solve the cubic equation 


(€*)PA, a7 (€*)PAy zB eA, + Xo —A = 0, (2.29) 


for No = 2777/4, Ag = Qn? 4, A, = = - 
to this same approximation is 


(WT) Mis = ar (Wy To. Vl 1 +4(W, 4 ml [(W, LO) mn ] 
_ a T,),,(1+7:15 x 10-74) (2.30) 


A, = 2:30 x 10-%. The heat transport 


i, 





m 


from (1.18), (2.18) and (2.19). In figure 2 we have plotted the total heat 
transport, divided by its value at Ay against A/Ay for the A, -approximation, 
for the A,-approximation from (2.27), and for the A,-approximation from 
(2.16). The two additional curves labelled «? and A, (ZANL) will be 
discussed shortly. 
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Effects of X, and zero-average non-linear terms 

Perhaps the most interesting aspect of the curve for A, is that it has 
brought the heat transport back very nearly to the A,-curve in the range 
Ayo <A < 3A). This alternating character of the various approximations 
to A makes the observed linear relation between heat transport and A more 
understandable. One can anticipate that the A,-curve, like the A,-curve, 
will diverge to the right at some A greater than 3A), while the A, 9-curve, 
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Figure 2. Heat transport vs Rayleigh number for the second, fourth and_ sixth 
approximations, for the fourth approximation with variable a®, and for the 
sixth approximation without zero-average non-linear terms. 


like the A,-curve, will be more nearly parallel to that for A,._ This alternation 
also suggests that a more satisfactory formulation of our problem should 
be sought permitting an expansion A = A(e) in some monotonic sequence. 

The percentage error of the A,-curve at various A can be estimated by 
comparing the magnitude of the various terms in the expansion of WT 
(or equivalently of 8/8,,). From (1.18), (1.13), and the preceding equations 
for W, and T,, 














H 


WT = &W, Ty + <2(We Tz + Wo Ty) + e4(Wo Ty + We T+ W, To) +.) 


{2 sin® wz + €?(1-098 x 10-*)sin zz sin 372 


ll 

m 
J 
— 
—_ 
— 
o™~ 
Ba 
~ 


+ €4[1-43 x 10-6 sin? 372 — 2-22 x 10> sin zz sin 372+ 


+ 2-38 x 10-5 sin wzsin 573}. (2.31) 
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‘Therefore 


+ (5:49 x 10-32 — 2:3 x 10-5€4)cos 42 + 1:26 x 10-*€4 cos 672 +...1. 


For A = 3A9, equation (2.29) has the root e? = 57:7, whence 

8/8, = 1+1-301[0-801 cos 27z + 0-160 cos 472 + 0-041 cos 67z+...], (2.33) 
suggesting that at this value of A the error due to neglected terms is less 
than 2°,.. This error is of course much smaller at smaller A, but will rise 
rapidly for A > 3Ap. 


Figure 3. Comparison of the mean temperature gradients predicted by the second 
and sixth approximations for free rolls at A = 3A . 


An interesting consequence of relation (2.32) is depicted in figure 3. 
The curves 8/8,, for A, from (2.17) and for A, from (2.33) at A = 3A, are 
drawn from a boundary to the middle of the fluid. Despite the negligible 
change in heat transport (given by the gradients at the boundary), the 
distortion of the finite amplitude disturbance has prevented the gradient 
from becoming negative. This is probably a property of the convection 
at all A, but no general proof that 8/8,, > 0 has yet been found. 

Another matter of some interest is the effect of the zero-average non- 
linear terms 4 and L(M) on the transport of heat. One can easily recompute 
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Ay, Ay, and Ag, neglecting # and L(M), and compare the resulting transport 
with that found in (2.29) and (2.30). From (2.13) and (2.24) it is seen 
that A, and A, are unchanged. However, a redetermination of W, and T, 
leads to Ag = 1-44 x 10-3 from (2.26). A new solution of the cubic (2.29) 
with this value of A, determines the heat transport curve labelled A,(Z4ANL) 
in figure 2. One concludes that, to this approximation, neglect of the 
zero-average non-linear terms increases the amplitude of the predicted 
heat transport. 


Amplitude effects on the horizontal wave-number 
Up to this point in the calculations we have considered only the growth 


5, its optimum value for a lowest Aj. 


process of the free roll when 2 
However, it is improbable that the distortions of the growing disturbance 
are not accompanied by a change in the basic horizontal wave-number. 
Indeed the expressions for Ag, A,, and Ay (see (2.11), (2.14) and (2.24)) 
have been written as explicit functions of z? so that we may investigate the 
effects of changing x on the heat transport. 

We have as yet no criterion to select the physically realized finite 
amplitude solution from the manifold of steady solutions permitted by the 
equations. In first-order theory the criterion is to find the lowest value 
of A,, for this disturbance is the ‘first’ to grow. But for A greater than 
this lowest Ay a variety of disturbances of different «* could grow, each 
corresponding to different A>Amin. In a following section we shall 
establish a ‘relative stability’ criterion to select the realized solution. 
Here, we shall investigate one extreme, namely, that solution which leads 
to the maximum heat transport. 

To the A,-approximation the heat transport is a maximum for minimum 
Ay, 1.e. at x? = 4, for all A. However from (1.19), (2.14) and (2.24), to the 
A,-approximation, 


— ( 2 _ 1/2) 
(WD), = 4 - 2 -[ (2) - oN) a (2.34) 
Ay Ay A, J 














where 
\ = 74(1+?)3 = 7?(1 +7)? a (1 +)? (32+ «?)? + 2(1+.«?)? 
rer 4” a 242 —— 4x2 (32 + a?)® — (1+ 2?)8 : 
We seek that #2 = 22(A) which leads to a maximum (W7),,.. From (2.34) 
this extreme relation is 
a OK 2K By y 
ee sani ae (2.35) 
oe a 20-2) —(WT)n, 
where K=A2/A,. An explicit solution for x? first will be sought for A 
close to (Aj)min Where z? will be close to }. We let 22 = }+A; then 
K = 6\,(1—aA + bA? +...) (2.36) 


and 








242 W.V.R. Malkus and G. Veronis 


where A, = (Ay)min, @ = 0-621 and 6 = 1-782 from (2.34). Hence 








(WT), = 6A,[1 — (1 —3x)!2+ AA — BA? +...], (2.37) 
where 1— x A—A, 
1 = a: ns =) >y c= 
oe: A 
Bn CTR TI See 
(i — gx)" 2 (i- jy" 
Then equation (2.35) for the extreme 2? yields 
A = A/ZB; (2.38) 
at A = 2A, for example, A = 0-0591 and 
(WT),,., 1-00672(W7),. 


Hence, for this extreme solution, z? does change appreciably even in the 
limited range of A in which the A,-approximation gives an adequate 
description of the motion. However, the change in heat transport is quite 
negligible and tells one that this change in x? plays an unimportant role 
in the energetics of our probem. 

‘Though certainly beyond the useful limit of the A,-approximation, the 
‘end point’ of the A,-curve has been computed as a function of x” to 


complete the extreme (W’7),,-curve labelled (x?) in figure 2. At this end 


point 

A= 4Ag/|Ag] +Ag= Ajnag? 
and 

(WT) mlaaa... = 4(Aimac — 0): (2.39) 

Hence, from (2.35) for the extreme (W7),, as a function of x, we have, 
using the values of A, and (W7),, in (2.39), 

CK / ex? 4 4 cr, Cr 0). (2.40) 
It is seen that this condition for maximum heat transport at the end point 
is identical with the condition for minimum A... from (2.39). A plot of 
Niue AGainst #, using Ap, Aj, and A, from (2.34), leads to the (A, H)-point 
given in figure 2. The value of x? at this minimum A... is 27/40. 


We have now progressed about as far as is practicable in the e-sequence 
for the simplest of the solutions to the Rayleigh problem. ‘The two important 
properties of the non-linear terms which control the amplitude of the 
steady disturbance have been exhibited. These are the distortion of the 
disturbance from its initial form and the distortion by the disturbance of 
the mean temperature field. ‘This latter plays the dominant role in 
amplitude determination for the free rolls. In the following section we 
shall establish the relative importance of these two ettects for the general 


rectangle and the hexagon. 


~ 


3. ANALYSIS FOR THE THREE-DIMENSIONAL CELLULAR MOTION 


‘he greater complexity of the general rectangle and hexagon plan-forms 


prohibits an analysis as extended as that for the rolls. However, the 
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non-vanishing of the first zero-average advection terms causes even the 


A,-approximation for these plan-forms to exhibit many of the properties 
found for the roll-cell in the A,-approximation. 


General rectangle 
We start this analysis with the general rectangular cell with two free 
boundaries. A complete solution to the first-order equations is 


W, = 2V2coszlx cos zmy sin 72, 


Ty) = No 2wv2 cos mlx cos zmy sin 72, 31 
Uy) = —22lx-* sin mlx cos mmy cos 72 ee 
Vy) = —2V2ma~ cos lx sin mmy cos 73, 
Ay = m(1 + a?)3/a?, No = 77(1 + «?)?/a, P+? = ee. 
From (1.5), ve and (3.1), 
Vilog = —87°N, a ?m*l*{cos 2zlx + cos 27my}sin 272, \ (3.2) 


Log = 873 /?m2x-4(1 + x?) [cos 2alx + cos 27my sin 272. z. J 
Hence A, = 0 as before and, from (1.20), 


168) = “e| . 1+ 3 OF cos ane 4 


° 


Ww 
— 


+41+ : Cs COS 2rmy |sin 2a (3. 
The solution of (3.3) is 
W, = (Cin, cos 27lx + C,,, cos 27my)sin 272, (3.4) 
where 
823x-2m2l? No{1+4o0-1(1+/)/(1+2?)} 
647°(1 + /*)3 — 42°PA, ‘ Cin 
From (1.21) and (3.4) 


T, = (D,,, cos 2alx + D,,,, cos 27my)sin 272, 


Cin = — (m, l) = Ci 





-_— 
w 
wn 

~~" 





where 
Ng Cin — 27-Ng mPa? 
Din = “FAC TP) : Din (m, 1) = Dyyr- 


From (1.9) and (1.10) expanded for U, and V,, 

e( M - oo ey = o(M, Jon Cx = 0, 
so that 
1 = —(C,,,/1)sin 2zlx cos 2zz, 


V,=-(C (3.6) 


) 
)/m)sin 2xmy cos 27z. } 


m 


To determine the first finite amplitude results for the general rectangle 
we must construct A, from (1.23). Here 


A; = i 0 Vi I olW, 0 To- W, Tolm i Wy Vi(hor +hyy)+ 
+a W, V?(Lo1 + L10)}m/ (Wo V2 Wo)ms (3-7) 
Q2 











jm 
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while from (3.1), (3.5) and (3.6) the lengthy forms of h;; and L,; lead to 


[Vin + Tao) 





a == (Di, + PD, .)= 4N,D(e,lm) (3.8) 


(W, V2 Won ial 
and 
2 = -_ 
[Wo Sa bon n+ Lw)In = SS (MPC, + PC) = 3No C(o, J, m). 
a V2 Wo) a4" 
‘Therefore 
A, = A,(0, m, 1) = $N,(1+D+C/o), (3.9) 


and from (1.18) and (1.19) this first etfect of « on heat transport is 


2(A —Aj) 








(Hi T),, = function of (o,/,m) = - (3.10) 


1+D+Cjo’ 
For the special case of a square cell /? = m? = $2?; from (3.4), (3.5) and 
(3.8) 





8x? 1 1 “ 
1+D+C/co=1 Ta8)(G2a8 = ne ee Ino)? (3.11) 
where x= (1+ $a?)/(1+2°) 
For the ‘limiting rectangle’ / > 0, m? > «?, we have 
14+D+C/o=1+4; (3.12) 
that is, W, = U, = V, = 0, from (3.4) and (3.6), but 
T, = —(N,/27)cos 2a/x sin 27 


from (3.5). 
Figure 4 is a plot of }(WT),,/(A—Ay) as given by equation (3.10), for 
values of //m between these two extremes, with «? = |, and for several 





values of co. ‘The dependence of (WT),, on x will be discussed shortly. 


Finite amplitude effects 

Perhaps the most important aspect of these finite amplitude effects is 
that the heat transport for the limiting rectangle does not coincide with that 
for the roll. Neither rolls nor long rectangles can satisfy the (distant) 
lateral boundary conditions in real convection. However, a long rectangle 
is certainly the more realistic limiting form to compare with squares and 
hexagons. 

A second aspect of this initial amplitude convection is that squares 
transport more heat than any other rectangle when o > 0-8*. Most gases 
have a value of o very close to this critical value and may exhibit a far less 
decided choice of plan-form than liquids. However there is some evidence 
that these strong effects of o on heat transport disappear at large Rayleigh 
numbers (Malkus 1954 a). 

As for the roll plan-form, one expects that the value of ” can change 
as the rectangle grows in amplitude. We investigate this effect for the 


* When o < 1 (cf. o 140 for mercury), the radius of convergence of the 
A = Ne, ao!) series which we have evolved is very small. In such a case it is better 
to non-dimensionalize V by writing V = 5V’ d and obtain the series \ = A(e, o). 
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case of the square and, again arbitrarily, seek that value of «? which leads 
to a maximum value of (WT),, atA > (Ay)min = A... From (3.10) and (3.11), 
epee 2 \ -1 
For values of A near Ay, « will be near }. We write A= «?—3}, whence 
Ay = A.(1 + 4A?2/3 +...), 
v = §(1—4A/15 + 8A?/45+...), 



































——- 2(A—X,) [1 —4A, A?/3(A—A,.) +... J x 
(WT), = —= od a (3.14) 
I (1—cA+dA*+...) 
Or 
oon hits. asacscah 
WT Mm Py = = 
SPRY OSE ~ 
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Figure 4. Heat transport for the general rectangle with a? = } at several values of o. 
where ; 100 ) : 
. F= 1+ a i+ See i aren ng 
473 \ 250 Ja~ 
100 = = me 9 
c= = . (0-7124 + (0-06450-1 + 0:2674c0-*), 
473F 
Oe care ie ) 
d = bi (07800 + 0-09276-1 + 0:29000-2). 
4735 
Therefore ¢(WT7),,/cA = 0, when the optimum value of A is 
A 4 A 1 
Aopt = —— es + —— (d-—c*)} 
m A = A a 3 A —A, 4 
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and [((WT)mlopt = 2(A—A,)(1 + dcAopt)/F. (3.15) 
For o = o, ¢ = 0-1245 and d = 0-1362. When A = 2A,, Aopt = 0:0428 and 
(WT), Jopt (WT),, 4?=4 = 1-00267. 


For o = 1, c = 0-1561 and d = 0-1589. When A = 2A., Aopt = 0-0532 and 
the above displayed ratio is 1-00415. Hence, for those values of o for which 
the square cell transports the most heat, this transport is only negligibly 
increased by an optimization of 2?. However this optimization causes 
a 10°) increase in x? for A = 2A 


o/\y- 





The optimum 2? equals } for the limiting rectangle since from (3.12) 
A, is independent of ?. 

The range of A in which these A,-approximations for the general 
rectangular cell are useful is apparently even more restricted than it was 
for rolls, i.e. A, < A < 2A,. This was established by the extremely tedious 
computation of A, for the square with x? =. It was found that here 
A, = 9-185, indicating that corrections due to A, become important for 
A—A, > 0°5A,. 


A, for hexagons 

Before attempting further interpretation of these results we will 
investigate the initial finite amplitude convection of the hexagonal 
plan-form. Past experimental literature gives one the impression that 
Rayleigh-like convection is realized only as hexagons. However these 
observations of plan-form were made with rigid bottom and free top 
boundaries. ‘The one experimental study of plan-form made with 
symmetric boundary conditions will be discussed shortly. 

The complete solution to the first-order problem for a hexagonal 
plan-form and two free boundaries was first given by Christopherson (1940), 
and is as follows 
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where L = 4/3 is the non-dimensional length of one side of the hexagon. 
Averages over the hexagon can be taken over the region 0 < x < v3L/2 
and 0 < y < L,2, which covers one-twelfth of the plan-form and is the 
smallest symmetric segment. 
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Possible variations of x? from its initial value x? = } have played such 

a small role in the initial heat transport of the general rectangle that we 

simplify the following analysis by choosing x? = } from the start. Then 
from (1.5), (1.8) and (3.16), 








Vi hoo » LN (db, + +h, »)sin 27z ea 
b (3.17) 
) ee — (32? 2)(4, T dy )sin 272, 
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= 5 eae ae pmo 
= + > ) 
d,= 2cos i cos — co —— 5 dy d, = 0, Vids - — (37/2), 
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Here, for the first time in our study, Ap, and Ly, have a horizontal dependence 
which is not orthogonal to the plan-form of W,. ‘This unusual property 
of the hexagonal solutions will be important in the case of one rigid and 
one free boundary, where due to the asymmetry the s-dependence of Mop 
and Ly) and the s-dependence of IV, are not orthogonal. However in this 


case of two free boundaries A, = 0, since (sin zz sin 27z),, = 0. ‘Therefore, 
from (1.20), 
f(W,) = $73 N,[(1 + 3/0), + (1 + 11/30), ]sin 27 (3.18) 
with the solution 
H ' (93 47)(( 19, + C ,» by )sin 273, (3 19) 
where 
e __i+3/o_ ins en ; 
729/4—Ac/7' - 1331/4—3Ac/z 


From (1.21) and (3.18), we have 


o «ffi. 3 es a SS - —— 
ie mal 3(F« a + (= C,- 5) 2 fpin2z. (3.20) 


The comparison of IW, and 7, for the hexagon with W, and 7, for the general 
rectangle assures one that the o-dependence of these initial finite amplitude 





distortions is similar. Therefore, a determination of the hexagon A, for 
a» x will permit an adequate comparison of hexagonal and general 
rectangle initial heat transports. Consequently the computation of A, 
from (3.7) is considerably simplified, for we need not compute Lo, and Lo. 
Then 


[Wo V2(hoy + Ayo) Im = — (222?/No)[To(hor + M10) In (3.21) 


from (3.16). Now (Th),, = () in general, and therefore 








(To hoa) = 9; (Ty oo + To(Ror +210) In = 9, ete; (3.22) 


hence 


[Wy V2(ho1 + hy) mn = (227?/No)(T; hoo) m- 
From (3.22), (3.21) and (3.18), 


7 73 2 22 0.7. 
(Tn n= — RCo + (ARO =) = — ntC. (3.23) 


we 
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Then from (3.7), (3.23) and (3.16) for ¢ = x, 


Ao = ANy—(T yhoo) = 4Nof1+4(C),. | = (972/4)(1-379). 
The initial finite amplitude heat transport for the hexagon from (1.8) 
(1.19), and A, above is (Hi Ti. 1-45(A—A,). In figure + we see that this 
heat transport is smaller than that due to the square and larger than that 
due to the limiting rectangle. 


Mixed horizontal plan-form 

‘There ts yet another class of solutions to the first-order equations which 
may grow to finite amplitude, viz. the class of all possible linear combinations 
of the individual plan-forms. We will consider only the single case of a 


square (S) and a general rectangle (G). ‘The first-order solution will be 


W, V2(d. 1 rd,.)(1 -r?)—-12 sin wz \ 2¢ sin 73, (3.24) 
where 

dy 2 cos(7ax y 2)cos(7%1 v2), d? = b 

¢,, = 2cos mlx cos zmy, + nt = a, 
and 0 <r < » is a free parameter determining the fractional contribution 
of ¢,. to the total plan-form. In the case z? = }, a= x, we can quickly 


determine A, by the method used for the hexagon (see equation (3.23)). 
From (3.24), 
Nog = 7No(1 +7?) Ago (S) + r7tog(G) 4 
t2r(l = ls m)cos( + l)ax cos( +m) sin 273}. (3.25) 
Hence 
W, = (1477) W,(S)+rPW,(G)+ 
+ (2rB/z3)(4+ A)~*cos(} +])zxcos(} +m)zysin2z3} (3.26) 
and 
T, = (147?)-47,(S)+7°7,(G) — (2rB/7)cos(3 + 1)zx cos(} + m)zy sin 273} 
where 
A=(3+l]P+(4imy, B= N,(1Flzm)(44+ A)?/[(4+ A) -AcA/z'}]. 
Then 
(Tha). : (1+r?) 21(T, hoo), (S) . (Ty Roo) n(G) — Lr2N, BIL F1F m)}. 


(3.27) 


Finally, from (3.24), (3.9) and (3.27) 

Ay = (14+ 7r?)-*fA,(S) + 277A,(1) +1A3(G)}, (3.28) 
where 

A.J) = 1Not+ i B(1F1F m). 

Hence A, has an extremum with respect to r at 

(7*)opt = [Ao(S)—A,(7)]/ [Ao(G) —A(1)], 
where 

(As)opt = [Aa(S)Ay(G) — A3(Z)]/ [An(S) + A(G) — 2A, (2)]. (3.29) 


This extreme is a maximum if 2A,(G)A,(/) > AS(/)+A3(G). Comparison 
of (3.9) for A,(G) with o = x, 2? = 3, and (3.28) for A,(/) establishes that 
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this inequality is true for all values of / = (}—m?)'*. We recall that an 
increase in A, decreases the heat transport. ‘Therefore, as one might have 
anticipated, this mixture of first-order plan-forms reduces the initial finite 
amplitude heat transport. 

Before a formal attempt is made to determine the physically realized 
solution from the vast degenerate set we have constructed, a relevant 
experimental study will be discussed. In 1935 Schmidt & Milverton 
made optical observations of the spacing of convecting cells between two 
rigid boundaries. Shining a collimated beam of light horizontally through 
water (o = 8) they detected a spatial oscillation in intensity due to the 
density variations in each cell. ‘The patterns obtained by rotating the 
beam through 180 in the horizontal contain one maximum in the spacing 
of the light fringes in the case of limiting rectangles, two maxima for squares, 
and three maxima for hexagons. At these maxima, from the initial 
a” = 3-13/7 (appropriate to the case of two rigid boundaries) the ratio of 
horizontal extent of each cell to the vertical spacing of the boundaries will 
be 2 for squares and limiting rectangles and 1-8 for hexagons. Schmidt 
& Milverton do not record any rotation of their beam to determine 
uniquely the plan-form. However, the average ratio of horizontal fringe 
spacing to cell height in their experiment was 2:1. They believed that they 
were observing squares. We can tentatively conclude that these cells 
were either squares or limiting rectangles, but not hexagons. 

The qualitative aspects of our heat transport computations for various 
plan-forms will most probably be preserved in going from the symmetric 
case of two free boundaries to the symmetric case of two rigid boundaries. 
For two free boundaries we found that the square transports more heat, 
the limiting rectangle less heat than the hexagon. Hence one anticipates 
that heat transport is a selective factor in establishing the observed flow 
from the many solutions. 


4. A ‘ RELATIVE STABILITY ’ CRITERION 

When seeking a realizable steady-state solution to the equations of 
motion one must verify not only that the solution formally satisfies the 
time-independent equations but also that it is stable against all infinitesimal 
disturbances which satisfy the boundary conditions. If it is not stable 
then surely there is at least one other solution to the equations. ‘This other 
solution may not be a steady one in the local sense, but for fixed boundary 
conditions is must be statistically steady; that is, integrals of moments of 
this solution over the entire field will be independent of time. 

The general problem of the stability of a known solution v and T can 
be immensely complicated. For example, years of labour have been devoted 
to simple two-dimensional shear flow {the Orr-Sommerfeld problem) and 
no exact solution has yet been found. Hence it is improbable that we could 
formally determine the stability of finite amplitude solutions even if we 
knew them exactly. 
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Fortunately the answer to a simpler question will suffice here. We 
will assume that we know the complete set. of steady and statistically steady 
solutions to the equations of motion. ‘Then any experimentally realized 
solution is contained in this set. We ask, are any of these solutions stable 
against those infinitesimal disturbances that have the form of the other 
solutions? If there is an affirmative answer for just one of these solutions, 
we will have determined the realized solution. If there are several solutions 
which are ‘relatively stable’, then either these solutions can be realized 
separately (i.e. they are metastable) or a broader class of infinitesimal 
disturbances must be considered to remove the remaining indeterminancy. 

In formulating such a question we have hoped that a detailed knowledge 
of the individual solutions and their interactions would not be required 
to answer it. Indeed we seek some simple integral property of the solutions 
which will single out the stable one. 

From (1.1), (1.2) and (1.3) in dimensional form, the equations of an 
arbitrary disturbance v’, 7” on the field v, T are 


(0/ot —KV?)T" v.VI’-v'.VT-v.VI, (4.1) 
(c/ct — vV?)v’ vy. Vv'—v’. Vv—v’'. Vv’ —- VP’ P,»tyl'k, (4.2) 
and V.v =0. (4.3) 


These are the general equations which determine the stability of v, T. 
We first construct the disturbance power integrals (see (1.6) and (1.7)), 


(7/2Bn)2(T’2) /Et = (y/Bm)[«(T’V2T"),, + (WTB), —(T’'V VT) yn] (4.4) 


1 A(v’.v’),,/et = v(v’. V2v’),, +7(W'T’),, —(v' ov’. Vv (4.5) 


where T has been replaced by T+ 7 and B 6T/dz. If v and T are 
stable to the disturbance v’, 7”, the right-hand sides of (4.4) and (4.5) 
will be negative. One may interpret (4.4) as the equation for the balance 
of entropy* associated with the disturbance 7’. Equation (4.5) is the 


m 


equation for the balance of kinetic energy of the disturbance v’.. The first 
term on the right-hand side of (4.4) is proportional to the loss of entropy 
by thermal diffusion; the second is proportional to the gain of entropy 
through interaction with the initial mean temperature gradient 8; and the 
third is proportional to the gain or loss of entropy through interaction with 
the initial field of temperature fluctuation. In (4.5) the first term is the 
loss of kinetic energy of the fluctuation v’ by viscosity, the second is the 
production of kinetic energy by convection and the last is the gain or loss 
of kinetic energv through interaction with the initial velocity field. 

We will now restrict the class of disturbances to be considered to those 
which have the form of steady or statistically steady solutions to the basic 


* In this irreversible process of steady convection we regard ‘ entropy ’ as propor- 
tional to and as a shorthand for the mean square of the local fluctuations of internal 
energy. We prescribe no (thermodynamic) role for this ‘ entropy’ which cannot be 
deduced from the basic disturbance equation (4.1). 
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equations; that is, 

v = Av, I" = BT, (4.6) 
where 4 and B are arbitrary amplitudes and v,, 7, are correct solutions. 
Then from (4.4), (4.5) and (4.6), 


Bt Uae (+7) 
where 





= —(T VT )n/(TP)m—— Le= (WTB) fr (TP 

f,= (v,T,.VT),,; I,= —[v(v,. V2v1) mn +hy)/(¥1- V1) mo 

p.=(.%. ¥¥) I,= y(W, T1)/(¥1-V1) 
The solution of (4.7) will decay with time if 
[411 + 13)? — HI — La Ly) ]'* -— 3 + Jy) < 0. (4.5) 


But (/,+/,;) contains the two large positive-definite diffusive dissipation 








m? m* 


terms and is invariably positive. ‘Therefore decay will occur if 
i,t, > [gly (4.9) 
Since v,, 7’; are solutions, 
K(7, VP; \mn 4 CW, l By) = 0, 
»(¥,- VV) m t¥(Wy Ty) =O. 


m 





(4.10) 


With the use of (4.10) the stability criterion (4.9) becomes 

[W, T1(81 — Bm tho tte Ti Bim/Y¥(M Tym > 9- (4.11) 
The basic equations for v,, 7; allow the triple product integrals f, and f,, 
to be expressed as double product integrals, that is, 


Ii = (v - 1 V2 )y2 =- (Tv. VT y)m | 
Soe ee I (4.12) 
f , = (v. dv, /at—vv. Vv, —-yWT),,. | 





All the double product terms of (4.12) vanish in horizontal integration 
if the field v,7 is developed from a fundamental periodic function 
orthogonal to the fundamental function of v,, 7;. This is the case for all 
different solutions which can be generated by the iterative method of the 
previous sections (for A < 24A,). Hence in our problem f,, and f,, vanish. 
To determine from (4.11) which of the many solutions is stable we must 
find that solution for which [W, T, T,(B;—B)], > 0, where the index 7 ranges 
over all solutions but the one considere H A physical interpretation of this 
conclusion is that the stable solution will produce more entropy per unit 
time from the mean temperature field of any other solution than it does 
from its own mean field. 

The stability criterion can be further simplified and given another 





«B+ VT for all solutions, then 


physical interpretation. — H= 


J 
[W;, Ti - Bn = BB;— BF) a 0. (4.13) 
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Condition (4.13) is not satisfied if (8*),, > (8?),,, for on applying Schwarz’s 
inequality we have (8?)*, > (f?),,(8?),, 2 (88;)%,. Therefore a necessary 
condition for stability is that 


(B*),, > (Bn (4.14) 


that is, the stable solution has a greater mean-square gradient than any 
other solution. If there is only one such solution, the criterion (4.14) 
uniquely resolves the formal degeneracy. If several solutions have the 
same maximum (§?),,, they are either metastable or a larger class of 
infinitesimal disturbances must be considered to decide which is realized. 

Yet another physical meaning can be given to the maximum (f?),, 
criterion. From (4.10), 


(4.15) 


m \ KD 


which states that for the stable solution the rate of dissipation of kinetic 
energy minus a quantity proportional to the rate of increase of entropy 
by thermal diffusion is a maximum. Equation (4.15) may also be written 


(62/82), -1 = NAF, 
where V=(WT),,|«B,». F= (WT?),, (WT) —1, (4.16) 


‘The quantity V depends only upon the amplitude of the convective heat 
transport. ‘The quantity F depends only upon the z form of the heat 
transport. But all the solutions for A, found in §3 have identical z forms 
for WT. Therefore the only stable solution found in $3 is the one of 
maximum convective heat transport. For o > 1 this solution is the square. 

A summary of our conclusions on relative stability is as follows. 

(a) Maximum heat transport is the stability criterion for all symmetric 
A, solutions. 

(b) Maximum (§?),, is the stability criterion for all solutions generated 
by the iteration of orthogonal linear forms (which includes all in this paper). 

(c) Equation (4.11) is the most general stability criterion for vertical 
convection. 

It is of some interest to relate the total potential energy of the convecting 





fluid to these stability conditions. In general the potential energy is V = gpz. 
Hence 
; go paz _ go a2 z= 00,, am 
(V), 2. | sp dz rie fa dz dz, (4.17) 
d , 2 d - di2 Joc 2 


where = is now measured from the mid-point of a symmetric convective 
field. From our density-temperature relation 
Cples’ = py a8, (4.18) 


so that 


(4.19) 
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The expression within curly brackets is zero when f = §,, and approaches 
unity as 8 approaches zero in all but the boundary regions. In steady 
convection 














so that | | 
— Vp) 4B (WT) 2  # aT 
i a Po ¢ m 1 = m. 1 ae Ks: TR Iz q . 4.20 
Pn 12 | 2B, | d)-asd (WT), | nee 


In comparing two solutions of identical form, i.e. with identical values of 
WT/(WT) 


will produce a minimum (V),,.. This conclusion is in keeping with the 


m equation (4.20) asserts that the one of maximum heat transport 
stability condition for A, deduced from (4.16). However, the general 
condition for relative stability is not the requirement that the total potential 
energy be a minimum though such a relation might have been supposed 
from considerations based on equilibrium mechanics. 


5. EXTENSION TO RIGID SURFACE BOUNDARY CONDITIONS 

In this final section we will discuss the determination of finite amplitude 
etfects for two rigid boundaries and for one rigid and one free boundary. 
We can hardly expect the rigid boundaries to modify the qualitative 
properties of the free boundary convection, for the important physical 
processes and the symmetry remain unaltered. We can expect quantitative 
changes due to the additional constraint on motion near the boundary. 
However, the asymmetry introduced by the rigid-free boundaries can 
modify the dependence of initial growth on A and, as experiments bear out, 
can lead to hexagons as the preferred plan-form. We shall find that the 
significant qualitative effects of these boundary conditions can be deduced 
without the final lengthy numerical computation for the A,. 

The complexity of these computations is perhaps the best measure of 
the limitations of the finite amplitude method studied in this paper. In 
conclusion we will turn from them to discuss alternative methods for 
determining finite amplitude effects. 


Rigid boundary analysis 
Proceeding just as in the simpler case, one writes the first-order solution 
for two rigid boundaries, as given by equations (2.9), (2.10) and (1.21), as 


‘ 
W, = > A; cosh 2,2 4(x, y)= We, | 
+ 5.) 
je (9.1) 
la al 6 _ 5) 9 
Ty) = a? > (4u2 —a?)A; cosh 2,2 4(x, y) = td, 
J 
where < is now measured from the mid-point, —} <2z< }, a@= 2?r?, 
A,|2 A, and py; are the complex numbers determined by Pellew & Southwell, 
A, is chosen to normalize W,, and ¢ is the normalized plan-form. In 


general 


— Vi hog = 2[ — (V2 g2)W at/ez + (— V2p)t OW) ea], (5.2) 
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Therefore 


Vr-h a : A, A; w,[(— Vid) (4u? — a*)? + (— Vid) (4u? — a?)?] x 


00 hon _ 


Lo =4> > A, A; u;{ (¢? — b*)(u? — vw?) + T?(4p? - a*)]|sinh 2(p;+ 1;)3, 


(5.4) 

where up; = — p_;, A;= A_;, Ay= 0. Hence, due to symmetry alone, A, = 0, 
because [cosh 2; z sinh 2(,;+;)s],, = 0. We let 

W, = 3) > P,,;sinh2(y;+ p,)z, (5.5) 


where — Vi P,;; = a? P,,;.. Then from (1.20) for W,, 
[(4(u; + w;)* —az}F4 az Xo | Pi; 
: 2A, 4; 41 (- Vi 6?) (42 —aq*)?+(- V2 Wh?) (4? — a)? + 
+ 4a?o "[ — Vi(* — fe?) — 4 (4, + w;)?(G? — be?) (uF — v7) 4 
+ ao [ — VET? — 4(u;+ w;)?T?] (4u?—a?)j,. (5.6) 
As we saw in $3, k sums over the three orthogonal functions cos 2z/x, 
cos2z7my and cos2z/xcos27my generated in hy and Ly, by the general 
rectangle, and over the two orthogonal functions ¢, and ¢, of (3.18) for 
the hexagon. ‘The 7, determined by this W, is written 
T, =1)> >G,,;sinh 2(u; + »;)2 (5.7) 
where 


a:G;.; P, A(t, + pL \- az\2 +2oG 14, A; ,\4(¢? - us?) (ps? p=) : 


Recall that 


dy = (Wy V2 Wo), [Wo To — (Wo Tom) Wo V3 Wo + Wo V2(hor + yo) + 
oT VE Ladin (5-8) 


‘me? 


01 
showing that the first finite amplitude results depend upon both the initial 
disturbance and the first distortion of this disturbance W,, 7,, U,, and V;. 
‘The only computationally difficult term in A, due to the (IW), T))-field 
alone is 

(WoT. Wo V2 Wo)m 


S A, A, A, A,(4u2- a2)? > 





| me + ty + pt, + fle) : sinh(p;+ 1; — pb, He) | (5.9) 
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which is the sum of 324 complex numbers. But the last two terms in A, 
caused by the distortion introduce 1944 additional complex numbers. 
‘The task of computing the numbers involved in (5.9) for several plan-forms 
can hardly be justified at this time. However, these numbers determine 
our only finite amplitude result rigorously applicable to a realizable 
experiment and could be used to check the physical validity of the formalism. 


Analysis for one rigid-one free boundary 
The case of one rigid and one free boundary leads to the first-order 
solution 
W, = > A;sinh 2p; z, (5.10) 
l 
where 0 <3 < } and A,/A, and yp; are known complex numbers as before. 
Comparing (5.11) and (5.1), we see that Vio) +o 1V*Lo, must be identical 
in form to (5.4), with the change that the index 7 runs from —3 to 3 while 
the index j runs from i to 3. Heyce 


(WoVE Wom As : [Wo(V? Roo ro IV? Loo) ] 





mt 


a? 3 A, 4, Ayn] 3 Prd | , 
r,j=1 h 


x [sinh 2p,2 sinh 2(u;+p;)2],, (5.11) 
where 
Pl, = {[(4u? — a2)*(— Ved) + (4p2—a2)2(— VEy2)] - 
— 4a?o" (uF — nF )[Vi(G? — f?) + 4(u; + By PP? — ¥")] - 


— @o (4? — a?) [Vi TP? +4(u;4+4;)?T7}},. (5.12) 


In contrast with previous computations of A,, the z-function average 


is not zero, while, as was noted in $3, P.;.¢ 


m 


ae eee, ee 
{sinh 2,2 sinh 2(p; 1;)=] 
will not vanish for hexagons. ‘Therefore a finite A, exists for hexagons 





with one rigid and one free boundary. Similarly (W, 7,+W, 7,),, will 


me 


not vanish and must be included in the heat transport. We will not proceed 
further with this particular evaluation of A, and (W, 7+ W, T7,),, but 
will enquire into the qualitative consequences of their existence. 

First-order theory is indeterminate with respect to the sign of the 
motion in the centre of hexagons. It could be either up or down. In 
contrast the sign of the amplitude of a rectangle is not a true degeneracy, 
since a change of sign leads to no observable change in the field of motion. 
The conclusion that A, is not zero uniquely removes this degeneracy in the 
hexagon, for A, determines the sign of the amplitude «. From (1.19), 
« = (A—A,)/A, for the initial convection, and hence if A, is positive € is 
positive and the motion is up in the hexagon centre, whereas if A, is negative 
« is negative and the motion is down. We have estimated A, by graphical 
methods for o = 8 and find that A, = 9-5. Hence the motion is up in the 
centre of hexagons as is observed. 
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One must now determine under what conditions the finite amplitude 
hexagon is preferred to the square or rectangle in the case of one rigid and 
one free boundary. ‘lo the A,-approximation, 

e*A, +€A, — (A—A,) = 0; 
hence 
e = (A,/2A5)[{1 + 4A,(A—A,)/A,}2? — 1), (5.13) 
where the second root is discarded since it does not vanish when A = Aj. 
The corresponding heat transport from (5.13) is 


(Wo Tom te(Wo Ti+ Wy T 5)n), (5.14) 





(WT),, = [(A—Ap)/Ao— €A,/A 


i L 


where we have kept all terms in (WT7),, necessary in the determination 
of A. Equation (5.14) asserts that the heat transport will be less than that 
due to A, alone for values of (A—A,) close to zero. However if 
\,(W, T,+ W, 7T,),,' (= QO say) is positive the heat transport will exceed 
the value appropriate to A, when 


ees 





S |Ay|(Wo T)1)2(Q/Ao)#? + (Wo Ty) Q. (5.15) 


2 m \X Y 
Since the value of A, for hexagons and the value of A, for squares will be 
comparable (certainly within a factor of two), then hexagons can become 
the preferred disturbance if A,(W,7T,+W,T7,),,! is positive and if the 
critical value of (A—A,) as given by (5.15) is small enough so that (5.14) 
is still valid. We have not yet determined the magnitude of (W,7, + W,T.),, 
but a graphical method of integration assures us that its sign is positive. 
Hence we conclude that hexagons will not appear as the initial instability 
with one rigid and one free boundary, but will appear for some finite (A—A,) 
with positive vertical velocity in their centres. There is some experimental 
evidence to support this conclusion. Benard (1901) describes the initial 
appearance of convection due to cooling a fluid from above as ‘‘a disordered 
cellular regime ’’ which becomes a steady field of hexagons after the cooling 
has continued a short time. 


Amplitude determination by integral techniques 

The lengthy task of computation needed to determine finite amplitude 
effects for realistic boundary conditions encouraged us to consider less 
exact methods. Stuart (1958) describes the use of power integrals 
to determine the amplitude of any disturbance whose form is assumed to 
be a good approximation to a correct solution. Stuart’s study was concerned 
with finite amplitude processes in shearing flow, but he has discussed with 
us the application of this method to thermal convection. We will briefly 
outline our impressions of the virtues and limitations of the thermal power 
integrals as tools for finite amplitude study. In particular, we will compute 
the heat transport for the case of two rigid boundaries which can be compared 
with the experimental results. 

In the steady state the kinetic power integral (1.7) is a homogeneous 
relation between the amplitudes and forms of v and 7. If v = Av’ and 
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T = BT’, where T” and the W’ component of v’ are normalized, then (1.7) 
in non-dimensional form may be written 

B/A = v(-v’. V?v’),,/y(W’T’) (5.16) 
On the other hand the thermal power integral equation (1.6) is an 
inhomogeneous relation since 8/8, from (1.13) also contains the amplitude. 
Equations (1.6) and (1.13) may be written 





me 


AB = (\—A,)/S, (5.17) 
where 
s.r, (E212. 252 
lel a 1 


Equation (5.17) for the convective heat transport is like the observed linear 
relation between heat transport and A, since S and A, will be constants for 
any particular form of v’ and 7”. If we use the form of the initial Rayleigh 
instability of the case of two free boundaries for v’ and 7”, then from (2.7), 


WR. 1, Sah Ady 


m 


and AB = 2(A—A,), (5.18) 


independent of the plan-form of the initial disturbance. Equation (5.18) 
is identical with our results for A, in the case of roll cells with free boundaries. 
This is understandable, for the free roll is the one case in which U,, V,, 
W,, 7, distortions and o effects did not influence initial growth. Since 
distortions and o effects uniformly tend to reduce amplitude we may hope 
that this use of the power integral sets an upper limit on heat transport. 
Other inhomogeneous integral relations may be deduced from the heat 
equation which lead to different but lower values of the heat transport 
at the same A and for the same choice of v’, 7”. For example we may write 
the non-dimensional heat equation from equations (1.5) and (1.13) as 
AW =h+[WT-(WT),,]W- VT L. (5.19) 
Then squaring both sides and averaging we obtain the inhomogeneous 
relation 
2 = (L2),,/(W?),,, (5.20) 
which must be satisfied by any correct solution to the problem. If we 
again use the form of the initial Rayleigh instability for the case of two 
free boundaries for v’ and 7”, but for simplicity use the roll plan-form 
(since (5.20) depends upon plan-form), then 


AB = (202—)2)!2—),. (5.21) 


0 
Thus we have derived a different law for heat transport from that given 
by the power integral. In general, equation (5.21) yields a lower heat 
transport for given A than does (5.18). Only for A = Ay do the two formulae 
agree since it is only at infinitesimal 4B that this choice of v’ and 7” is a 
correct solution. In spite of this difficulty, the power-integral heat transports 
seem to bear too remarkable a resemblance to the observations to be 
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fortuitous, as we will now show by a computation of the Rayleigh function 
for rigid boundaries paralleling (5.18). 

Pellew & Southwell have proposed an approximate solution to the 
first-order characteristic functions (for rigid boundaries) which is easier 
to work with than the correct solution, 

W”’ = N!2{sin naz + Pcosha(z— $)+ On(z— })sinha(z—$)]¢, (5.22) 
and 
T’ = v2¢sin 72, 
where ¢ is the normalized plan-form, N?? is a normalizing constant, 
P = (-4921 and QO = 0-3416 are constants chosen to satisfy boundary 
conditions. Using this form and relation (1.21) to fix B/A, equation (5.17) 
leads to A, = 1713°7, S = 1/1:51, and 
AB = 1:51(A— 1713-7). (5.23) 
‘The correct value of A, for the infinitesimal disturbance is 1707-8, indicating 
that (5.22) is an adequate approximation to the correct infinitesimal solution. 
The heat transport law (5.23) is in such good agreement with the observations 
reported in table 1 that one is tempted to extend this use of the power integrals 
beyond the range of steady convection. As discussed in the Introduction, 
a second instability occurs at A = 18000, leading to aperiodic convection 
and an abrupt change to a new linear law of heat transport. In all, five or 
six transitions appear in the data up to A = 10*, where ‘fully turbulent’ 











| | | 
| | | | | 
vssiatil A A | sp | ABIA = uf | 
Fansition | Theoretical | Experimental | Theoretical |Theoretical mee | 
mental 
l 1 708 1 700 +4 80 | 1-51 0 0 | 
2 17 600 | 18 000+ 1 000 | 1:72 1-37 1:3+0-1 | 
3 61000 | 55 000+ 4500 | ieRs | 2563 2:3+0-1 | 
4 | 170 000 170 000 + 15 000 | 1-9 4-29 3:4+0°1 | 
5 | 411 000 425 000 + 20 000 | 1-9 6°87 4-5+0:1 
6 | 855 000 | 860 000 + 30 000 | 2-0 | 8-52 5-7+0-1 | 
| | 





Table 1. Convective heat transport at the experimental transitions computed from 
equation (5.24). 

conditions seem finally to be achieved. If, as suggested by the success of 
the power integral, self-distortion of a disturbance is not very important, 
then perhaps the interaction between several ditferent disturbances plays 
a small role in the determination of heat transport. ‘To test this possibility 
we can determine the total heat transport which would result from the 
independent growth of first, second and higher Rayleigh-like modes of 
instability. ‘That is, 


AB = ¥ (—Ac)/S, (5.24) 


l 
where the S; determine the slope of the growth curve for the mode 7 and 
As, is the value of A for instability for that mode. Equation (5.24) has the 








= = 
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form of the observed heat transport up to A = 10%. The values of the S, 
and A, have been computed in the same way as S and A, of (5.23) with the 
approximate functions of Pellew & Southwell for each instability. The 
results are tabulated in table 1 and compared with observations in water 
made by Malkus (1954a). In dimensional form AB/A is the ratio of the 
convective heat transport to the heat transport in the absence of convection. 
For large values of A, S; |! approaches its free surface value and has been 
set equal to this value for A greater than 800000. Equation (5.24) agrees 
quite well with experiment at the end of the steady convective range, but 
by A = 55000 at the end of the first aperiodic range there is an error of 
approximately 15°... This percentage error continues to increase with A 
through all observable transitions. Equation (5.24) cannot be correct in 
the range of beyond observable transitions since it leads to an incorrect 
gross dependence of heat transport on A (giving AB/A ~ A!‘ instead of the 
observed variation as A!3). However it is rather remarkable that (5.24) 
works as well as it does. This suggests to us the possibility that aperiodicity 
permits the several disturbances to be more or less independent which 
certainly could not be the case if they were steady superimposed motions. 
The quasi-independence due to the aperiodicity could permit a greater 
release of potential energy and hence lead to a more stable field of motion 
than any other. Such thoughts are of course only speculation and a formal 
theory for the onset of thermal turbulence has yet to be found. 


6. CONCLUDING REMARKS 


To summarise the conclusions: the initial finite amplitude convection 
is determined primarily by the distortion of the mean temperature field, 
secondarily by the self-distortion of the disturbance; the number of formal 
steady-state solutions increases with increase of Rayleigh number beyond 
the critical value; however this degeneracy is removed if there is only one 
finite amplitude solution which has a greater mean-square temperature 
gradient than any other solution. More particular conclusions are that: 
square plan-forms are preferred to hexagonal plan-forms in ordinary fluids 
with symmetric boundary conditions; the initial stages of convection are 
markedly altered in a fluid of small Prandtl number; the effect of distortion 
of the disturbance on the thermal field prevents the mean gradient of the 
temperature from changing sign in the mid-regions of the fluid; the zero- 
average non-linear terms reduce the amplitude of the disturbance; the pre- 
ferred horizontal wave-number increases with increasing Rayleigh number. 

The approach to finite amplitude processes given above is applicable 
to those convection and shear flow studies with adequately solved stability 
problems. ‘The method can determine those conditions which permit a 
steady finite amplitude disturbance to occur before the appearance of the 
infinitesimal disturbance. Hence it is a tool to extend the study by Meksyn 
& Stuart (1952) of the onset of disturbances in channel shear flow. The 
method is of particular value in the resolution of explicit or hidden 
degeneracies in a stability theory. Hidden degeneracies exist when there 
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are velocity-dependent forces other than advection which play no role in 
the criterion for instability. ‘This occurs for example in the study of 
geomagnetism. ‘here one wishes to determine when an electrically 
conductive fluid heated from below will select a joint magneto-convective 
state of motion rather than a pure convective state. 

However, beyond the initial finite amplitude effects the technique we 
have described here can become prohibitively tedious. In addition the 
e-sequence will fail to describe the preferred field of motion after a second 
instability has occurred, though this sequence may still converge to a set 
of formal steady-state solutions. We continue to search for an adequate 
descriptive framework which can include the second instability and in 
particular can deal with aperiodicity. 

This study bears a relation to two previous works on fully developed 
turbulence. A paper on turbulent convection (Malkus 1954b) and a 
paper on turbulent shear flow (Malkus 1956) were based on the assumption 
that the most stable field of flow would be that one which released the 
maximum amount of potential energy per unit time. ‘The section of this 
paper on relative stability suggests that a maximum mean-square temperature 
gradient (equation (4.14)) is the more appropriate stability criterion. The 
only qualitative prediction of the earlier works which appears to be altered 
by this change is the structure of the boundary regions and this causes a 
reduction in the turbulent transports. These modifications will be 


discussed in a following paper. 


It was our good fortune that Dr J. T. Stuart was visiting the United 
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SUMMARY 

‘This paper deals with the stability of a two-dimensional laminar 
jet against the infinitesimal antisymmetric disturbance. ‘The curve 
of the neutral stability in the (x, R)-plane («, the wave-number ; 
R, Reynolds number) is calculated using two different methods 
for the different parts of the curve; the solution is developed in 
powers of (zR)~! for obtaining the upper branch of the curve 
and in powers of «R for the lower branch. 

The asymptotic behaviour of these branches is that for 
branch I, x—2, c--2 for R-»«%; and for branch II, 
R ~ 1:12%-!?, ¢ ~ 1:20%2 for «+0. Some discussion is given 
on the validity of the basic assumption of the stability theory in 
relation to the numerical result obtained here. 


1. INTRODUCTION 

A large part of the theory of hydrodynamical stability, as it has been 
developed in the past, is concerned with the question whether an essentially 
parallel flow is stable or not against infinitesimal disturbances. Within 
this limitation the theory has been successful in predicting the stability 
characteristics of a number of typical laminar flows, its results having 
been supported satisfactorily by experiments. 

The application of this theory, however, has hitherto been confined to 
flows with at least one solid boundary, and little is known about the stability 
of boundary-free flows such as a laminar jet, wake, and mixing region between 
parallel flows. ‘This may not be surprising in view of the fact that the free 
flow, having at least one inflexion point in its velocity profile, is supposed to 
be highly unstable so that its critical Reynolds number is very small. For 
such a small Reynolds number the usual methods of asymptotic approxi- 
mation for large Reynolds number become ineffective and a new method 
of approximation must be introduced. The more fundamental difficulty 
of the boundary-free case, however, lies in that the basic assumption of the 
hydrodynamical stability theory that the undisturbed flow is approximately 
parallel is not satisfied for smaller Reynolds numbers, because then we can 
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neglect neither the lateral velocity component nor the change in the 
longitudinal velocity profile. ‘Thus, the investigation of the stability of 
free flows under the framework of the existing stability theory appears to 
lead to self-contradictory results. ‘The present calculation on the two- 
dimensional laminar jet seems to confirm this anticipation; the critical 
Reynolds number is found to be only 4-0, associated with the wave-number 
a = 0-2. 

The stability of the laminar jet has been already dealt with by Curle 
(1957), using a new method of approximation which was devised by McKoen 
(1957) and himself, especially for investigating the stability of free flows. 
‘Their method is essentially to reduce the fourth-order disturbance equation 
to a second-order one, using some physical considerations. However, to 
neglect the fourth-order term in the equation does not seem allowable for 
small wave-number even when the Reynolds number is fairly small. More- 
over, Curle’s way of approximating the unknown stream function by the 
linear combination of two inviscid solutions is quite arbitrary and a different 
choice of the combination factor may lead to a different value of the critical 
Reynolds number. 

The principle of the calculation adopted in this paper is quite simple. 
Solutions are expanded into inverse power series of «R for large value of «R, 
and in ascending series of zR for small «R. In practice, solutions are 
calculated up to the second- and third-order term for the above series 
respectively. Unfortunately, owing to the slow convergency of the series, 
we have not been able to work out the entire curve of neutral stability in 
the (x, R)-plane. But the branch of the curve for x0 is sufficient for 
determining the critical Reynolds number. 


2. FORMULATION OF THE PROBLEM 


The basic flow under consideration is the two-dimensional laminar jet 
ejected from an infinite line orifice. ‘The velocity distribution of this flow 
was calculated by Schlichting (1933) and Bickley (1937), giving the result 


u= Uysech'(), 
: 2.1) 
v= (#R,)-*U je sech>{ 2 ) — tanh{ 2 | 
ie "Ld L L}f’ | 
3M2\"8 M \-18 
mits = . oem = ) ai 


x 
and Mx\!? 
R, = (=) : (2.3) 


M= | u® dy = const. (2.4) 


— 2 


where 





represent the Reynolds number and the kinetic momentum flux of the flow, 


respectively. 
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It may be seen easily from (2.1) that for large Reynolds number R, 
the ratio v/u (= O(R, *4)) is small, so that the flow becomes nearly parallel. 
Moreover, since from (2.1) (eu/dx)/u = O(x-"), the variation of u with x 
is also small for large values of x which are associated in general with large 
R,. ‘Thus, the fundamental premise of the hydrodynamical stability theory 
that 

du/dx = 0, v = 0, (2.5) 


is approximately satisfied in this flow for large R,. Now let us assume 
that the condition (2.5) is satisfied, expecting R,. in the present calculation 
to be fairly large. It will be found later that the final result of the calculation 
does not always confirm this assumption, but this will be discussed in § 7. 
We take U, and L defined by (2.2) as the representative velocity and 
length of the flow, respectively, and make all variables non-dimensional 
with respect to these characteristic quantities. Denoting the dimensionless 
coordinates again by (x, ), the non-dimensional velocity profile is given 

from (2.1) by 
U(y) = sech? y. (2.6) 


It is convenient to define another Reynolds number by 


R=.U,L/» (2.7) 
which is related with R,. by 
R= (5) 3Re. (2.8) 
The relation (2.8) shows that the discussion of the previous paragraph on 
the condition (2.5) is equally valid when R., there is replaced by R. 
By virtue of Squire’s theorem (Squire 1933) we need consider only 
the two-dimensional disturbance, whose stream function may be defined by 


b= dene", (2.9) 


where x, being real and positive, represents the wave-number of the 
disturbance, and c is in general a complex constant. ‘The real part c, of ¢ 
represents the phase velocity of the disturbance whereas xc,, c; being the 
imaginary part, is the amplification factor. According as the values of c¢; are 
positive, zero, and negative, the basic flow becomes unstable, neutrally 
stable, and stable, respectively. 

Substituting (2.9) into the Navier-Stokes equations and neglecting the 
non-linear terms with respect to ¢, we obtain the Orr-Sommerfeld equation 
for the disturbance : 


(U—c)(¢" —2°6)— U"b = — (¢" — 2276" + 26), (2.10) 
iaR 
where the primes represent differentiation with respect to y. 
The boundary condition for the disturbance is given by the requirement 
that all disturbance velocities must vanish at infinity, that 1s, 


d' = 4b = 0) aty = +. (2.11) 


When U(v) is an even function of y, as in this case, we can consider the 
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antisymmetric and symmetric disturbances separately. Since the anti- 
symmetric disturbance is, on physical grounds, supposed to be more unstable 
than the symmetric one, and this is the case for flows with solid boundaries, 
we shall deal here only with the former disturbance. For antisymmetric 
disturbances, the boundary condition (2.11) reduces to 


d'(ax) a.h( mo) = (0) : 6” (0) = (), (2.12) 


and we need consider only the semi-infinite range y > OU. 


3. SOLUTIONS FOR LARGE aR 
The analytical property of the solutions of the Orr-Sommerfeld 
equation (2.10) for large values of zR has been investigated in detail by 
previous researchers (sce Lin 1955, ch. 8). Its four particular solutions 
are classified into two inviscid solutions and two wvtscous solutions. While 
the former two have non-vanishing values throughout the flow field for the 
limit «R—- ~, the latter two vanish almost everywhere except inthe immediate 
vicinity of the solid walls. We need not therefore consider these viscous 
solutions in the present problem, where there is no solid boundary, so long 
as we are dealing with large values of «R. Foote & Lin (1950) clarified this 
point mathematically and concluded that the effect of viscosity enters the 
stability problem of the unlimited flow only through the higher approxi- 
mations of the inviscid solutions, and that the boundary condition at the 
origin ¢'(0) = 46”"(0) =0 (for the even solution) reduces to #’(0) = 0. 
We therefore deal here only with the inviscid solutions which are obtained 
in the form 
= 8) 
dv) = > (a2R)*d™(y; «, €). (3.1) 
Substituting (3.1) into (2.10) and equating terms of the same powers of «R 


on both sides, we obtain the following equations : 


(U c)(d 92h) as U’" dO ~ 0, | 
(l c)(d 2h) — U"g™ (fp —Y — 2a2hit -! + tgin—)), 


for n> i. 

The solution of (3.2) is first calculated in three separate regions of the 
tlow field: (i) the region far outside of the core of jet, (i1) the region near 
the critical layer at which y = v,, and (iii) the region around the centre. 
Then, on joining these solutions analytically at two boundary points of 
these three regions, we obtain the solution over the whole flow field. 

(1) Outer solution. Yor large y the velocity profile may be written in 
the torm 

. 
U(y) = 4e-2" > (-1)™(m+ 1)e-™, (3.3) 


so that 


—— = > A, e-*, (3.4) 


m= 1 
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where 
16 


Cc 


Ax (1—2c)A,, 


al-> 


A, = 1 (16 —40c¢+27c?)A,, Ay= + (16 —48¢ + 46c? — 16c*)A,, 
a 


9 


Ss 


Substituting (3.4) into (3.2) and taking account of the first two boundary 
conditions of (2.12), we obtain the following solutions for large y: 
) 
a0 | 
gf ~ Do — Cie ay 2 Ane amy } 
lies (3.5) 
«a 
fp) ~ OW = oe enw = b,,, en7am ve | 
m=1 J 


where C; is an integration constant and 


A, A,a,+Ae 

















oor l, = ’ i <r ee wae, 
a ay 4(1 +a) ay 8(24 xs 
A,a,+A,a,+ Az 
a, 
‘ 12(3 +) 
5, = Sa + %)ay, 
ys 
64i +a)? ) 
b, . | 1, 5, oa (2+«)?a, o ay |, 
~ 8(2+% l c 
| 16t {. .2 x 16 - 
4, = i234 5 | Ais + Ayb, + = ~~ + a)ay 4 = (2+ «)*a, + 
‘. a”) 
; 3 (2 —~c)(1+«) a} | 
/ 


(it) Critical-layer solution. In the neighbourhood of the critical layer 
y=y,, Where U(y,) =, the solution is expanded in a power series of 
(v—y,). Substituting the Taylor series expansion 


i | ere 2 
U-—c= U,(y-y,)+ 71 UN(Y¥—Ye)* + + 


into the first equation of (3.2), we obtain a pair of independent particular 
solutions for ¢ as follows: 


00 
pf?) pa yo = c. > Sy —y,)", 
m=1 


a 





(3.6) 
— ut... | 
pd ~ FO = af 2 Lmn(¥—Ve)" + TU! I log(y —y, Mr» 
where C, and C, are integration constants and 
: U? : a2 = ye Uiv 
= I, o= =, i ae ae eis oll Ae pee A . 
h Ss 2uU, Is 6 6U fs 7) Ss 24U, 
9, = 1, oy = (), pe eae Ue = (Z). 
ss " — 2 wae U' 
eU" UU! Lae - ae 3 
7 18 U; 6U?? Zu, 8\ Ut] 


The way of choosing the proper branch of log(y—y,) in ‘IQ’ has been 
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established by earlier workers (see Lin 1955, p. 130) with the result that 
if L’’ < 0, as in this case, 

log(v—v,) = log y—y, for y—y, > 0, 

log(y—y,) = log'y—y, +ix for y—y, < 0. 


The second-order solution 4" is obtained from the second equation 
of (3.2) by quadrature and two particular solutions are expressed in terms 
of the 6s as follows: 


Pw PO = PO | MPD dy—YO | M,'V¥® dy, (3.7) 
where : ” : 
M, = —i(U—c)3(FO" —2aFO + o2PO), j= 1,2 
The determination of the proper logarithmic branch must be carried 
out in the same way as for 4“. It should be noted here that the solution 
‘VD is O\(y—y,) 7} in the neighbourhood of the critical point, so that the 
l 


) loses its meaning there. Since, however, we are not concerned 


solution ‘I 
with the eigenfunction itself but only with the eigenvalues of x, R and c, 
a good approximation at the joining points of the solutions (y, and y, in § 4) 
is sufficient for the present purpose. In fact this is attained in the present 
case, because the join ng points are fairly distant from the critical point 
and the coefficients of the singular terms in 4“)) are small for small (c—c,)/c, 
(c, is the limiting value of c for zR—> «), and therefore the singular terms 
are insignificant at the joining points. 

(iii) Inner solution. For small value of y the velocity distribution may 
be written in the form 


U-c= > F,,y™, (3.8) 
where 


Fy =1-c, F,, = (—1)m22%m+b(2%+0 _ 1) 


2m +1 
2 : 
(2m +2)! In bb mM 2 l, 
and the B’s are the Bernoulli numbers. 
Substituting (3.8) into (3.2) and taking the third boundary condition 
of (2.12) into account, we obtain the solution in the neighbourhood of the 


origin as follows: 


d ws OO). c. > h,, y™, 





m=0 : (3.9) 
dP~QV= Cy > ky", | 
m—1 
where 
l toa l ee ; shaken 
hee I, h, ap (Fo + 2F,), h, = Tar, (eh + 12F,)+ 02F, Ay}, 
h., = a (2? Fy + 30F 3) + (2? F, + 10F,)A, + (2? Fy — 10F;)hy} 
1 


hy = oe ((27F3 + 56F,) + (a? F, + 28F3)h, + 0? F hy + (2?Fy — 28F,)h3}, «5 
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; 








k= - oF (24h, — 47h, + «'), 

— — { — i(360h, — 2402h, + ath,) + 02F, h,}, 

— = {_ i(1680h, — 60x2h, + athy) + (a2F, + 10F,)k, + (a2Fy— 10F, ks}, 
= 0 


4. BOUNDARY VALUE PROBLEM FOR LARGE aR 
We now proceed to join the outer, critical-layer, and inner solutions 
analytically at some appropriate points, y, and y,, say, each located on the 
boundaries of the three regions. Noting that the system of equations (3.2) 
is of the second order, we can easily see, extending Foote & Lin’s (1950) 
analysis, that the effect of the viscous solution does not enter the eigenvalue 
problem and that the analytical continuation of the solution is simply carried 
out by making the functional value and the first-order derivative of the 
solution continuous at y, and y, instead of doing so up to the third-order 
derivatives. If this process is carried out the solution thus obtained will 
satisfy all the boundary conditions of (2.12), because the outer and inner 
solutions derived in §3 already satisfy the respective boundary conditions 
for y-> x andy = 0. 
First, joining the outer and critical-layer solutions at y, (¥, <¥, < &), 
we have 
@'(y,) — Wir) + (Cs/C2) V2) 
O(y,) WiC) + (C3/C2)¥ 2091) - 
‘The similar relationship between the critical-layer and inner solutions 
is obtained by the continuation at y, (0 < y2 < y,): 
Q"(y2) _ ¥i(V2) + (C3/C3)'¥2(v2) 
Q(y2) Vs (y2) + (C3/Ce)¥o(¥2) 
Eliminating C;/C, from (4.1) and (4.2), we obtain the eigenvalue equation 
relating the parameters «, R and c: 
E(a, R,c) = 0. (4.3) 
We are particularly interested in finding the condition for the neutral 
disturbance. Putting c; = 0 in (4.3) and solving the complex equation, 
we can express « and R as functions of c 
a= a(c), R = R(c). (4.4) 
The limiting case of the problem for «R— «% was solved by Savic 
(1941). According to his results, 
a= 2, U"(y,) = 0, 





(4.1) 


(4.2) 





so that 
c= U(y,) = 3, 
and the solution (for the antisymmetric disturbance) is given by 
¢@ = const. x sech? y. 
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lor large «R the corresponding value of c may be expected to be not 

far from the limiting value c, = 3. ‘Then, for a number of the prescribed 
values of c in this atime we obtain an equation connecting x and R, 

and we can find the roots of x and R by trial and error. During this process 
we expand the coefficients of the critical-layer solutions into powers of 
c—c, and neglect the terms of higher order than (c—c,)°. It is found that 


s 





























l 
| c x | R | 
| 
§ (¢ 's) 2 ( xg) | wD 
| 0-66 1-97 137 | 
| 0-65 1-93 | 270 | 
| 
0-64 1-88 161 
0-63 1-84 116 | 
(0-62 1-80 91-9 
0-61 1:76 76:7 
0-60 1:71 66°5 | 
| 
‘Table 1. 
2 
Branch I 
(a6 
1 
Present result 
/ 
/ — Curle’s result 
Branch II ~- 
| 
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Figure 1. The curve of neutral stability in the (a, R)-plane. 


the real solutions are obtainable only for those values of ¢ smaller than c, 
and that the approximation of the present solution becomes unsatisfactory 
for c < 0-6. In the process of continuation of the solutions, the points 
v, = ¥,.+0-2752 and vo = 0-5503 are used as the joining points, and the 
solutions are calculated up to the fourth term for ®, the fifth term for Q, 
and the seventh term for ‘Y and ‘Q respectively. The second-order 
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solutions are calculated to the same order of approximation as the corres- 
ponding first-order solutions. The error caused by neglecting the 
remaining terms of each series is estimated to be at most 5°,,.. The numerical 
results are tabulated in table 1 and shown graphically in figure 1 (branch I) 
and figure 2. 
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Figure 2. The curve of neutral stability in the (c, ~)-plane. 











5. SOLUTIONS FOR SMALL aR 
We now proceed to calculate the solutions for small values of «R. 
Equation (2.10) may be written as 
(D? — «?)(D? — «? + iaRce)d = 1aR{U(d" — a?) — U"d}, (5.1) 
where D = d/dy. In the region far from the core of jet, that is, for large 
value of y, the velocity profile U() almost vanishes as well as all its spatial 
derivatives. ‘Thus, the right-hand side of (5.1) being negligible in this 
region, the four particular solutions of (5.1) reduce to 


wn 


dh = ety, exsy, where f? = «?—iaRce. (5.2) 


An approximate solution may be obtained by solving (5.1) with one 
of the solutions in (5.2) substituted on the right-hand side. Better solutions 
are derived by repeating this iteration process and we can express a formal 
solution for ¢ in the form 

© 
d(y) = > (iaR)"d™(y; a, B), (5.3) 


n=O 


where the ¢’s are related to each other by the following equations: 


D2 — x2)(D? — B2)46 = 0, ) 
ae ee | | (5.4) 
(DP? — a*)(D* — B*)g = U(g*-) — ag") — Ug", asl. | 

The uniform convergence of the solution (5.3) for ali «R can be shown 
as follows. Equation (5.1) together with the boundary conditions (2.12) 
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For large ~R the corresponding value of c may be expected to be not 
far from the limiting value c, = 3. ‘Then, for a number of the prescribed 
values of c in this neighbourhood we obtain an equation connecting « and R, 
and we can find the roots of x and R by trial and error. During this process 
we expand the coefficients of the critical-layer solutions into powers of 
c—c, and neglect the terms of higher order than (c—c,)*. It is found that 





























| | | 
| c x | R | 
| | | 
| ic) | 2 (a,) | L 
0-66 1-97 | 737 
0-65 193 | = 270 
0-64 1-88 161 
} 0-63 1-84 116 | 
| 0-62 1-80 | 91-9 
0-61 1-76 76-7 H 
0-60 1-71 | 66°5 | 
= 
Table 1. 
2 — | 
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ox 
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Present result 
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Figure 1. The curve of neutral stability in the (a, R)-plane. 


the real solutions are obtainable only for those values of c smaller than c, 
and that the approximation of the present solution becomes unsatisfactory 
for c < 0-6. In the process of continuation of the solutions, the points 
vy, = ¥,+0°2752 and y. = 0-5503 are used as the joining points, and the 
solutions are calculated up to the fourth term for ®, the fifth term for Q, 
and the seventh term for ‘Y® and ‘VQ? respectively. The second-order 
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solutions are calculated to the same order of approximation as the corres- 
ponding first-order solutions. The error caused by neglecting the 
remaining terms of each series is estimated to be at most 5°,,.. The numerical 
results are tabulated in table 1 and shown graphically in figure 1 (branch I) 
and figure 2. 
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Figure 2. The curve of neutral stability in the (c, «)-plane. 











5. SOLUTIONS FOR SMALL aR 


We now proceed to calculate the solutions for small values of «R. 
Equation (2.10) may be written as 


(D? — ?)(D? — 2? + ine) = ixR{U(¢" —02$)—U"d}, (5-1) 


where D = d/dy. In the region far from the core of jet, that is, for large 
value of y, the velocity profile U(y) almost vanishes as well as all its spatial 
derivatives. Thus, the right-hand side of (5.1) being negligible in this 
region, the four particular solutions of (5.1) reduce to 


h = e**Y, etsy, where B? = «?—iaRc. (5.2) 


An approximate solution may be obtained by solving (5.1) with one 
of the solutions in (5.2) substituted on the right-hand side. Better solutions 
are derived by repeating this iteration process and we can express a formal 
solution for ¢ in the form 


d(y) = > (taR)"d™(y; 2, B), (5.3) 
where the ¢’s are related to each other by the following equations: 
D2 2 D2 2\4(0) — 0), ] 
( x”) (D? — B®) L 64) 


(D? — «?)(D? — B2)h™ as U(d” Wy” _ gy 2din 1) — U" dh -1) n> E. 


The uniform convergence of the solution (5.3) for all «R can be shown 
as tollows. Equation (5.1) together with the boundary conditions (2.12) 
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is equivalent to the following integral equation: 
(Vy) 


| cy 
Ae-* + Beau + =| —e 2" | 
: 


(a ee 
U'e"d dy emt U'e-*"4h dy : 
o . x 
y "7 B? - 42 
ter? | (: + - 


i a 
- UL’ \e&"h dy + 
2p : 


‘y Jf) (Be OF 

+ eB | (« - 5 U, Je Bud, dy | 
= Ae~*" + Be? 4 

tia2R 5 E a (+ 


B)(m —y)} sinh} }(x—B)(m—y)! 
(a+ B) 








where A and B 
to 2 


}(%—B) 
+U(n) sinh} B(y — 


Y); = 
8 = H] b(n) dy, (5.5) 
are integration constants. 
é / (0) < ~ « 


Changing the variable from y 
2 1), equation (5.5) takes the following form 
d(z) = f(z)+A | K(z,6)d(C) dZ, 
where 


(5.6) 
ixR, f(z) = Az*+ Be*, 


cat ee (£) (: l 
ae tte z sia 
‘The iteration process suggested in the previous paragraph leads to the 
solution in powers of A 
d(s) =f(s)+A | K(s, 0 f(f) d+ > AM) K(z,C 
0 Raa 


2, fi)dl, | K(L1 2) dl... 
Rb hele 2) 
Now f(z) < |A\+|Bi, 
and since U(¢) = 4¢7(1 + C7)? 
U(¢)/¢| < 44, U'(¢)/¢ 
Then, if « < 2, A(B) < 2 


K(z,¢)| < v(2) : for 0 < C< 
where p - 


1 and N is a finite constant. 


v 





Hence 
- ; t N 
(" K(z, eye) dt |< GALAIBDA 
0 1 p 
so that 
> " , ” Sox, — : Al 1 R\ s 
| | Kis, 0,):@&, 0. | 'K(¢ > Sn f(%,) a, a1) (N 
~ O l—p m 
- (|A|+ Bi) N 
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‘Thus, the series (5.7) as “ as (5.3) ne uniformly for all values of 
\ (= iaR), provided « < 2 and A(f) < 
The system of equations (5.4) can lie, oer successively by the method 
of variation of parameters, leading to the result as follows: 
dp”) = ent, p) = ex! pf) baa “By fp) = eB: 
1 u 


gp) = Rl ay | Uer(d\" lt’ 4 ag\" ») dy 1 


(5.8) 


= ad 





9 9 
a2 + B? 
— eBy J ‘yy 8 wes > t pi - dy — 
\ 


2 Q2 | 
oo 22 gon) ] 


forn>1,j;=1,2,3,4. ‘The general solution of (5.1) is given by 
b = Cy $+ Crh. + C3h3 + Cy dy, (5.9) 
where the C’s are arbitrary constants. 
6. BOUNDARY VALUE PROBLEM FOR SMALL «R 
According to the first two conditions of (2.10), d, and ¢, are rejected 
and the lower limit of the integrals in (5.8) must be «. ‘Then, in order 
that the second two conditions of (2.10) be satisfied by the solution (5.9) 
with not identically vanishing C, and Cy, ¢, and ¢3 must satisfy the 
following condition: 
4:0) -4,(0) 
6/0) (0) 
which gives the eigenvalue equation between the parameters x, R and c. 
Substituting (5.8) into (6.1), the characteristic equation becomes 


| a+ B24 S (iaRYT(S,) — ¥ (iaRY"T(}g) | 


| n=1 n=1 


|= (0, (6.1) 


. ’ (6.2) 
S (iaRY"J(4) x2— BB+ ¥ (iaRY"IO(bs) 
1 n=1 
where ; 
I(g;) =2 | Ufsinh(xy)d\" + acosh(zy)h"~)} dy, 
i + (6.3) 


J(;) =2 i: U <sinh(By)d\"~? + ~ at 2 cosh(By)d\ "f dy, 
We work out the solutions up to the second approximation al neglect 
terms of ixzR of higher order than the third in equation (6.2). ‘Then we have 
— (x? — B?)? + in R(a? — B?){I(p,) —J(hg)} + 
+ (aR)PL(22— B)LI™(4,) F(a) 
+ [9(h,)I(b3) —1(¢ O(b,)] =0, (6.4) 


where the /’s and ./’s are calculated using (5.8) and (6.3). 
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In order to work out the integrations in (6.4), it is convenient to approxi- 
mate the velocity profile U(v) (= 4e-?”/(1 + e-*”)?) by a polynomial in e~*” ; 


m 


Uly) = > Gye". 
k=1 


In practice we adopt the following quartic 

U(y) = 3-87e-?4 — 6-50e—4" + 5-61e—*Y — 1-98e-*", (6.5) 
which fits the exact velocity profile quite satisfactorily as shown in table 2. 
Although the profile given by (6.5) does not satisfy the condition U’(0) = 0, 
this does not matter because all terms in (6.4) are expressed as integrals 
involving U alone. The complex equation (6.4) substituted from (6.5) 

















y | Exact U(y) | Approx. U(y) | 

@ és peti ee | 

0 1 1 

0-10 | 0-9972 1-0086 | 

()-22 0-9879 0-9973 | 

0-36 0-9689 0-9748 | 

0-51 0-9375 0:9372 | 

0-69 0-8889 0-8875 

0-92 0-8163 0-8163 

1-20 } 0-701 | O-7115 | 

1:61 | 0:5556 0:5557 

2:31 0:3306 0-3274 

| ee | 
| | Ud ; 3 09975 
| a | U dy 0:-6931 | 0-6861 

Table 2. 


is rather cumbersome and therefore we again expand its second and third 
terms in powers of B—«z (= a0, say) and neglect terms of higher order 


than O(c). ‘Then, equation (6.4) takes the form 

14+7R V(a, c)— R? W(2, 0) = 0, (6.6) 
from which R is given by 
W.. (6.7) 
where the suffixes r and 7 refer to the real and imaginary parts of the quantity, 
respectively. Eliminating R between (6.6) and (6.7), we obtain the equation 


en Ne 
ao os = if) 6.8 
1- EV, (x) WwW. =0, (6.8) 


from which o (or 8) is determined as a function of x. Substituting o(x) 
into (6.7), R(x) is obtained, and finally c(«) is calculated from the relation 


of 
[= - \] l + 3? 6.9 
iR ( said 
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A numerical calculation has been carried out for « = 0-1, 0:2, 0-3. For 
a > 0-4 the approximate equation (6.6) does not give a real root, and more 
laborious manipulation of (6.4) is necessary for obtaining the neutral 
curve for larger x. Numerical results are given in table 3 and figure 1 
(branch II) and figure 2. 











| x R € 

| O-1 | 4-21 0-018 
| (2 4-02 0-061 
| 0-3 4-39 0-121 
| 





Table 3. 


- 


7. DIscUssION 

It may be seen from figure 1 that for large Reynolds number R the 
wave-number « decreases monotonically with R, although the change in « 
is not remarkable for R > 200. ‘The same tendency of the neutral curve 
was noticed by Lessen (1950) in the case of the mixing region. On the 
other hand, the upper branch of all existing neutral curves for the bounded 
flows behaves in a quite different manner; « increases with decreasing R. 
It may therefore be inferred that the viscosity acts on the unlimited flows 
as a decaying factor, whereas it has the dual effect of decaying and amplifying 
the disturbances of bounded flows, although the two examples only may 
not be sufficient for so concluding. 

The lower branch of the neutral curve gives the critical Reynolds number 


R=4-0 at «=0-2, (7.1) 

which, according to (2.8), corresponds to 
R,, = 3-77. (7.2) 
The numerical values of the critical Reynolds number given by (7.1) and 


(7.2) may not be very realistic, because for such small Reynolds numbers 


like these the basic assumption of this calculation (2.5) is not satisfied very 
well. For instance, the ratio of the lateral and longitudinal velocity 
components of the main flow becomes, from equation (2.1), 

v/u = 0-5 x (non-dimensional function) 


for R, = 3-77, so that the second assumption of (2.5) can no longer be 
valid under such circumstances. On the other hand, the first assumption 
of (2.5) is not directly affected by the numerical value of the Reynolds 
number because the ratio (eu/ex)/u depends only on x. Rigorously speaking, 
therefore, all that we can say is that (7.1) gives the critical Reynolds number 
of an artificial flow with the velocity profile (2.6) which must be maintained 
by applying some sort of body force. ‘To treat the stability problem of the 
non-parallel flow in a more satisfactory manner seems to be beyond the 
scope of the existing theory of hydrodynamical stability. Having ourselves 


F.M. Ss 
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no new ideas at the moment for improving the existing theory so as to 
be able to deal with the more general case, we content ourselves with taking 
the values given by (7.1) and (7.2) as a qualitative estimate of the critical 
Reynolds number of the laminar jet flow. 

‘The asymptotic behaviour of the lower branch can be obtained exactly 
without using the approximate formula (6.5). Retaining only the lowest 
order terms with respect to « in (6.4), we obtain 


B(x+8)+itxzRix2—B+3B(x+ B)log 2! + a?R%4a— 28 —3(«— f)log2' = 0, 
(7.3) 


where use is made of the relations | Udy=1, | Uydy=log2. If 


we solve this complex equation with the relation 8? = 2? —1%Rc, we obtain 
the asymptotic behaviour of the branch 


a= is. Ba bie, 


(7.4) 
and c = {4(aR?)—! —2}a? = 1-20x?. 


With this relation (7.4) between «, R and c, it is confirmed that the higher 
order terms in 1xR which are neglected in (6.4) do not produce any term 
of the same order as those in (7.3), and therefore (7.4) gives the exact 
asymptotic behaviour of the neutral curve for x > 0. The asymptotic 
branch (7.4) is shown graphically in figures 1 and 2. 

The relation (7.4) shows that for x — 0, «R tends to zero whereas R 
itself increases indefinitely. This confirms the self-consistency of the 
expansion (5.3) of the solution in ascending powers of «R. It should be 
noted here that the so-called inviscid solution, which corresponds to the 
limit «Rk —- «x, has nothing to do with the lower branch of the present case 
which leads to the limiting ‘inviscid’ disturbance in the sense that R > x, 
xR —-0 tor x-»0. ‘This asymptotic behaviour clearly disproves Curle’s 
interpolation formula which expresses the unknown stream function as 
a linear combination of two inviscid solutions. It may also be interesting 
to note here that Pai (1951) found that the flow is unstable for the combination 
of large value of «R and small « He concluded from this that there is no 
lower branch of the neutral curve or the lower branch is very close to the 
R-axis in the («, R)-plane. According to the present result, only the latter 
of his conclusions can be correct. 

In figures 1 and 2 the neutral curve calculated by Curle is shown for 
comparison. His curve agrees fairly well with the present result so far as the 
upper region of the (z, R)-plane is concerned. ‘This is not very unexpected 
because his first assumption of neglecting ¢'° might be permissible for 
values of x which are not small and the possible error involved in his 


interpolation formula may also not be serious for small values of (c—c,)/c,. 
However, it may be too early to decide from this single example whether 
this agreement is entirely accidental or whether it gives some justification 
for McKoen and Curle’s new technique, at least, for the neglect of the 
¢'* term on the upper branch of the neutral curve. 
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SUMMARY 

The equations of motion for a chemically reacting gas in the 
absence of viscosity and heat conduction are set up. It is shown 
that the characteristic speed defined by this set of equations is the 
high-frequency limit of the phase velocity of sound waves as long 
as the reaction rate is finite. At infinite reaction rate (chemical 
equilibrium) the characteristics suddenly change to the low- 
frequency sound speed. ‘The nature of this transition is discussed 
in connection with a recent paper of Resler (1957). 


1. INTRODUCTION 


From the point of view of chemical thermodynamics there are no 
essential difficulties in writing down the equations governing the motion 
of a gas which is not in chemical equilibrium. ‘This is done, for example, 
in the theory of irreversible thermodynamics and in papers on laminar 
Hame speed. In these cases the main accent is on the chemistry. 

If one is interested mainly in the dynamics of the gas flow some changes 
in the choice of variables prove to be advantageous. A special case, flow of 
an ideal gas with lagging vibrational or rotational heat capacity, has been 
discussed from this gasdynamical point of view in an earlier paper (Broer 
1950). One of the results of this work was that the characteristics of this 
type of flow are determined by the high-frequency limit of the speed of 
sound. 

Apparently this fact has been sometimes thought to be slightly puzzling. 
In a recent paper Resler (1957) has given another treatment from which 
it would appear that this property does not hold. Asin the author’s opinion 
this treatment is not entirely satisfactory, it is proposed to return to this 
question, generalizing some of the work of our first paper. We shall pay 
especial attention to the relations with sound dispersion and characteristic 
theories, and since we are solely concerned with the interaction of reaction 
and flow, we shall neglect the effects of viscosity, heat conduction, and 


diffusion. Moreover, for the sake of simplicity, a gas in which only one 
reaction occurs will be considered. 
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2. EQUATIONS OF MOTION 
Under our assumptions the general equations of gasdynamics take the 
following familiar forms. ‘The continuity equation is 


e +t pdivv = 0, (1) 
the Eulerian equation of motion is 
Dv 1 
Dy’ r= (2) 
and the energy equation is 
= (Jo? + U)+ ~ div(pv) = 0, (3) 


where U is the energy of the gas per unit mass and the other symbols have 
their usual meanings. Introducing the enthalpy per unit mass H = U + p/p 
and using (1) and (2) we derive two equations equivalent to (3) 

a l 

= ly? 4 H) ey! 
Dr pot 
and DH 1Dp 


— (4) 


| 

| 

| 
On 


In steady flow equation (4) yields the equivalent of Bernoulli’s law 
$v? + H = const. on a streamline, but (5) is the most convenient form for 
most of our work. 

In the gasdynamical treatment it 1s advisable to consider H as a function 
of p and p. Moreover, H will depend on a parameter q which determines 
the chemical composition of the gas. ‘The energy equation (5) then can be 


H a + (H ;) - H ~*~. 0. (6) 
CC a ” Dt 
If there were several reactions we would have to consider a set of parameters, 
but we confine ourselves here to the case of one reaction. In order to have 
a complete set of equations, however, we still need an equation for the rate 
of change of g. It will be assumed that for each value of p and p, there 
exists one equilibrium value of g, which is denoted by 4(p,p). The 
tendency of the reaction will then be to make g tend towards the momentary 


written 





value of g. We now suppose that g can be chosen in such a way that 
Dq _ 
a 
is a good approximation to the rate equation. ‘lhe reaction rate z (a positive 
quantity) depends on the state of the gas and is therefore a function of p, p 


~4(q—4) (7) 


and, possibly, g. 
Once the functions H(p,p,q), 9(p,p) and «(p,p,qg) are determined, the 
set of equations (1), (2), (6) and (7) is sufficient to describe the flow. 
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3. DISPERSION OF SOUND 
As in the theory of sound propagation, we put u = eu,, Pp = py + €p, etc., 
substitute in the equations of motion and retain only terms of the first order 
in «. In order to obtain plane progressive waves we take all first-order 
quantities proportional to exp?(wt—kx). ‘The resulting set of homogeneous 
equations for the amplitudes can be solved only when its determinant 
vanishes. ‘This condition then yields the required relation between w and k. 
Remembering that gy = ¢(Ppo, po) since the undisturbed gas is in equili- 
brium, the indicated programme can be carried through by means of simple 
and straightforward calculations. ‘lhe resulting expression for the square 
of the phase-velocity of sound is 
w —H,(tw+%)—2H,q, 
Rk? (H,—1/p)(tw+%)+2H,@,, 
In this dispersion equation the suffix zero in the coefficients has been 
dropped. We can rearrange the equation by introducing limiting speeds 
for high and low frequencies 
ie -H i —(H,+H,4q,) 


and ¢ = ~ : (9) 
(H,,+H,4,)-—\/p 





(8) 





From the thermodynamical relation 


(=) (=*) 2a (=) 
Cp}p\op/s p op . 
it is seen that ¢ corresponds to the usual value for the sound speed when q 
is kept fixed. This is to be expected since the reaction will be frozen at 
sufficiently high frequencies. On the other hand, é can be obtained by 
supposing that equilibrium is maintained throughout, i.e. q = g. Writing 
H{p, p,4(p.p)} = A( p,q), we have 
f= 8,+84, H,, = H,,+ Hq, 

which proves this contention. 

Substituting (9) in (8), we obtain 


w c — Ce 

— = ?— —___, 10 

k2 1 +iwr ( ) 
where _ ae H, 








ac? H,+ Hq, 
At intermediate frequencies the phase velocity is complex, which means 
that there is absorption of sound. 

We note here for further reference that owing to the absorption the 
variable parts of pressure and density are not in phase. ‘The continuity 
equation (1), after linearization, reduces to 

lwp, — potkv, = 0, 


and the Eulerian equation (2) to 


: i 
lwv,— —tkp, = 0. 
Po 


Therefore p,/p, = w?/k?, which is complex. 
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In the usual theory of sound dispersion due to lagging heat capacities 
H=C,T+C°%, 
where C,, is the non-lagging specific heat and @ the temperature of the 
lagging degrees of freedom. If we now take 6=4, H = (C, R)(p/p) + Cg, 
it can be verified that this theory is a special case of the formalism of this 
section. 


4. CHARACTERISTICS 

The set of equations (1), (3), (6) and (7) is quasi-linear, which means 
that it is linear in the derivatives of the dependent variables 7, p, p and q 
but not in these variables themselves. A characteristic equation of this 
type involves only one specific combination of derivatives. We consider 
in this section only one-dimensional unsteady flow. A characteristic 
equation then can be interpreted as a law of propagation. Speed and 
decrement of this characteristic propagation depend only on the variables, 
not on their derivatives. 

In our problem, (6) and (7) are already characteristic equations as they 
stand. ‘The speeds involved are equal to v, and the equations therefore 
determine relations between rates of change of some quantities along the 
path of a fluid element. 

Using (7) and (9), we write (6) in the form 


oe. & eT q)- (11) 


Dt Dt 4H, 


‘The remaining two characteristic equations are linear combinations 
of the continuity and Eulerian equations. When (11) is used to write 
the former equation in terms of derivatives of p, it is easily seen that the 
required combinations are 


D..p + pe(D.,,v) = —ac?(H,/H,)(q—4), (12) 
where oe D - a 
Dt Ox 


This equation corresponds to characteristic speeds v+c. According 
to the results of §3, this implies that the characteristics correspond to 
propagation of a disturbance of infinite frequency. 

In the general theory of characteristics, which was not needed here 
since the required linear combination could be found by inspection, the 
characteristic speeds are found from a determinant equation. It is easily 
verified that this determinant has exactly the same structure as that required 
in the dispersion theory, apart from the terms in x, which do not occur 
in the characteristic determinant. ‘The algebraic reason for our result 
is therefore obvious. ‘The physical interpretation, however, has caused 


some difficulties. 
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5. MOTION NEAR EQUILIBRIUM STATE 
‘Lhe trouble with the interpretation of the characteristic speed lies in 
the motion near equilibrium. In this case, the change in condition of a 
fluid element is slow compared with a, and q—@ therefore remains very 
small. It is to be noted that this does not mean small Mach number or 
Mach number change. ‘The restriction is on the acceleration, not on the 
velocity. 
When q-~@ is exactly zero, the energy equation takes the form 
Aa Dp A Dp Dp 
v a v ag nas 
Dt Dt Dt 


or 2 Dp Dp | 0) 
Dt Dt ; 
Using this equation in the continuity equation, we obtain the linear 
combinations 
D,-p + péD.;7 = (), (13) 


which are again characteristic equations, but now with v + é as characteristic 
speeds. 

This result is not unexpected as the gas in equilibrium can be considered 
as a non-reacting gas with appropriate equation of state, since its condition 
depends only on the two variables p and p. ‘he question, however, is 
what happens when the reaction rate, although large, is finite. One could 
conjecture that the characteristic speeds would in some way change 
continuously from v7 +c to v+é, but this appears not to be true. As long 
as there is a finite reaction rate, the theory of §4 applies and 7 +c is the 
characteristic speed, although it would be expected that the motion tends 
to a solution of (13) with characteristic speeds wv + é. 

Before investigating this behaviour, we first point out that for motion 
near equilibrium q— 4, but not «(q— 4), is a small quantity since this motion 
is realized only when « is large enough. We have in fact, from (7), 


ee 


—a(q—q) = ne i Vp Dt 10 Di’ 
Therefore, there is at any rate no direct contradiction between (12) and (13), 
since the right-hand side of (12) does not tend to zero with increasing 
reaction rate. 
We shall now derive equations equivalent to (12) but with their left- 
hand sides identical with (13). We start from the energy equation (6) 
and substitute 


ee Be, eR ER 


Di Dt Dt! 2) = % De *% De * Det? 


In this way we obtain 


= = — — (gq q). (14) 
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When (14) is substituted into the continuity equation, we can form 
the linear combinations 


D.-p+pcéD.; eee (q—@). (1: 


Ji 
— 


‘These equations tend directly to (13) when the reaction rate x goes to 
infinity, because the right-hand side depends on (q¢—g) but does not, as 
in (12), contain a factor « 

In a numerical calculation procedure use of (15), instead of (12), could 
be advantageous when the reaction rate is sufficiently high. Nevertheless 
it must be kept in mind that the set (15) is not a set of characteristic equations 
in the usual sense. The reason is that the coefficient at the right-hand 
side does not depend on the variables only but contains derivatives. 
A consequence of this fact is that (15) does not have all the properties of 
characteristic equations. For example, discontinuities in derivatives do 
not propagate with the velocity +é but with +c relative to the gas, no matter 
what the reaction rate is as long as it is finite. Only when g = q exactly, 
do the derivatives drop out from the coefficients in (15) which then goes 
over into (13) and so assumes all the properties of a characteristic set of 
equations. 

The remark on the propagation of discontinuities of derivatives furnishes 
another clue to the situation. Let us consider a small disturbance, confined 
in a bounded region at some moment and having discontinuities in some 
derivative at the boundary of this region. When we linearize the equations 
and take Fourier transforms the spectra will extend to infinity. It is 
therefore not unexpected that the boundary of the disturbance will propagate 
with the high frequency sound speed c. 


AN ALTERNATIVE FORM OF THE EQUATIONS 
A very simple variant on the equations (12) and (15) has been proposed 
by Resler (1957). It is easy to show, starting from the one-dimensional 
continuity and Euler equations, that 


D.,p+paD,,,0 = 0, (16) 
where , _ Dp Dp Re 
ee te ie (17) 


‘These equations look like characteristic equations. However, it is not 
possible to express a in terms of v and the variables of state without using 
derivatives. One could eliminate, for example, Dp/Dt by means of (11) 
or (14) but then Dp/Dt will still persist. The equations (16) do not 
therefore have the properties of characteristic equations. ‘They might, 
however, on occasion be useful for numerical work, just like (15). 

It could easily be assumed that a would lie between c and é, but this 
would be fallacious; it is not even true that a always exists. ‘This would 
require that Dp/Dt and Dp/Dt were of the same sign everywhere, and we 
saw in $3 that in the simple example of a plane travelling sound wave there 
is a phase difference between the variable parts of p and p owing to the 
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reaction rate. ‘his entails a phase difference between dp/ot and cp/dt, 
which therefore change sign at different moments. 

In view of the engaging simplicity of (16) this result is somewhat 
annoying. It has, however, a sound physical basis. ‘The equations (12) 
essentially describe propagation of waves with velocity c and a certain 
attenuation owing to dissipation of mechanical energy which is given by 
the right-hand side. Equation (17), in connection with (16), means that 
one looks for the speed of propagation a of a point which has a constant 
wave amplitude as the propagation equations have no right-hand sides 
in this case. ‘This is possible only as long as the chosen value of the 
amplitude is less than that at the summit of the wave. 


7. DIsstPATION 
In terms of the usual thermodynamic variables, we can write 
dH T dS + dp p+O dq (18) 
and consider H as a function of S, p and gq. O = (cH/cq),,. is a quantity 
of the nature of a chemical potential. (Q is different from H, = (¢H/eq),,., 
as defined in § 2.) 

As the enthalpy H tends to a minimum for constant values of p and S 
we can calculate g from the equation O(g,p, S) = 0. Therefore, when q is 
close to g, O will be approximately equal to B(qg—q), where B = 0?H/cq* 
is a positive quantity since H must be a minimum for g = @. 


When, for given values of p and S, only one equilibrium value of q is of 


interest, O does not change sign at other values of g. Hence, O(q— 4) is 
always positive or zero. We now combine (5) and (18), obtaining 

DS De 

oe 4ahee 

Dt ~ Dt 

‘This equation corresponds to DS/Dt = 0 in ordinary gasdynamics, and 
it could be used instead of the energy equation. Substituting the rate 
equation (7) we find the dissipation equation 

DS _ 


i < «aie «th. ( 
! > xO(q-4) (19) 


It is seen, therefore, that the dissipation is always non-negative. 
When « is large, (19) is approximately 

DS _ 

—* 

which is then of order x~'. This yields the expected result that the 

dissipation vanishes in the limit of fast reaction. 


a i 





1 xB(q—9)?, 


‘This paper is Communication no. 89 from the Laboratory for Aero- and 
Hydro-dynamics of the ‘Vechnical University at Delft. 
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SUMMARY 


A new water tunnel, incorporating a slotted wall working 
section, was found to suffer from severe vibration. A theoretical 
explanation is given for this, together with experimental evidence 
gleaned from this water tunnel and a small wind tunnel. It is 
shown that the oscillations are hydrodynamic in origin and are 
associated with the slotted wall design. Consideration is given 
to methods of elimination or reduction of the oscillations. 


1. INTRODUCTION 


The slotted wall tunnel is well known as a tool for studying transonic 
flow. ‘The tunnel! choking and interference effects associated with a closed 
tunnel severely restrict the size of model which can be tested at transonic 
speeds; on the other hand, open jet tunnels are unsteady and the jet breaks 
up. ‘The slotted wall working section is intended as a compromise, the solid 
part of the boundary steadying the flow and the open part of the boundary 
preventing choking and reducing interference effects. 

Even in very subsonic conditions, these advantages are still valid. 
In a closed tunnel it is not desirable to use a model whose diameter exceeds 
about one-sixth the working section diameter because of interference effects, 
and an open jet tunnel is unsuitable if the length of working section is 
greatly to exceed the diameter of the jet. ‘The replacement of a closed 
working section by a slotted wall section either enables larger models to 
be tested (up to about one-third the diameter of the jet) or, at the design 
stage, permits the use of a smaller tunnel with consequent saving in expense, 
space and power. However, one present disadvantage of the slotted wall 
tunnel is the lack of design experience, so that there is an element of risk 
in embarking on such a design. Until this state of affairs is remedied, it 
is inevitable that the potential advantages of the slotted wall will remain 
unrealized. 

When new hydroballistic facilities were planned for the Admiralty 
Research Laboratory, a water tunnel was proposed to permit the testing 
of torpedo-like bodies of up to 10 in. in diameter. Since this would entail, 
for a closed tunnel, a working section diameter of 60 in., it was decided to 
take the risk and to design, instead, a slotted wall working section of diameter 
30 in., the length of working section being 15 ft., or 6 diameters. A general 
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view of the tunnel is given in figure 1; the total length of the upper limb, 
which includes the working section, is 71 ft. and would be almost doubled 
with a 60 in. diameter jet. A full description is given by Lever et al. (1957) 
from which source figure 1 has been adapted. Experience with the tunnel 
has given much design information, which, it is hoped, will put the design 
of a slotted wall tunnel on as firm a basis as the design of a closed or open 
jet tunnel. 

















Figure 1. General view of 30 in. water tunnel. 


Although the essential requirement is merely that the jet boundary be 
partly open and partly closed, attention has in the past usually been focused 
on walls slotted in the direction of flow. A sketch of the essentials of such 
a working section is given as figure 2. ‘The obvious relevant parameters 
are the number of slots and the ratio of slot width to slot spacing. ‘The 
designer, in addition to determining appropriate values of these parameters, 
also has the task of successfully guiding the jet into the diffuser, the criteria 
to be satisfied being a constant static pressure distribution along the axis 
of the working section and sufficiently small wall interference effects on the 
largest sized model to be tested. A certain amount of work ts available 
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on these subjects, both theoretical and experimental (e.g. Wright & Ward 
(1948); Davis & Moore (1953); Vandrey (1955)), and a small-scale 
mock-up investigation is appropriate. As a result of such an investigation, 
the slotted wall parameters and the diffuser dimensions were determined 
for the A.R.L. tunnel. At the time, it was expected that a satisfactory 
tunnel would result, particularly since the principles were first applied to 
a 12 in. diameter tunnel with apparent success. However, from the first, 
large tunnel vibrations were observed, the entire upper limb being set 
into motion. A certain reduction was obtained by anchoring the outer 
shell of the working section but the pressure fluctuations were still excessive 
both from the point of view of safety and of flow conditions in the working 
section. 








\ = f 

/ . T — ee! 

{ ae } mal al Ba . 

L aa ———L OS See Be DIFFUSER 
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Figure 2. Schematic diagram of a slotted wall working section. 


A return was made to the small wind tunnel, which is both more easily 
instrumented and more quickly modified during test. As a result, the 
basic explanation of the phenomena occurring was arrived at, but it was 
found, on investigating the water tunnel further, that certain differences 
in detail existed between the two tunnels. ‘This paper describes the basic 
phenomenon, which is the existence of a feed-back loop whereby particular 
disturbances may be much amplified, draws on experimental evidence 
from both tunnels to demonstrate this phenomenon, and, finally, describes 
briefly some possible remedies which have suggested themselves as a result 
of this work. 


2. OUTLINE OF MECHANISM OF INSTABILITY 


Self-sustained oscillations in a closed circuit demand that the total 
gain round the loop be unity and that the phase change round the loop be 
an integral number of cycles. In practice, the gain exceeds unity for small 
disturbances, and the amplitude of the final disturbance is determined by 
non-linearities in the system which ensure that gain is a decreasing function 
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of amplitude. Since losses exist in the circuit, an essential for self-sustained 
oscillations is that one element of the circuit should be an amplifier, i.e. it 
should be capable of drawing on a source of energy (in this case, the steady 
tunnel flow) and converting it to oscillatory energy. ‘This function is 
performed in the tunnel by the jet itself, bounded by the slotted wall. ‘The 
remainder of the circuit, the return path, controls the frequency of the 
oscillation by virtue of the phase change it introduces. ‘l'wo effective return 
paths have been found in our investigations; it is considered unlikely that 
others exist as their attenuation would almost certainly exceed the gain 
of the amplifier element. 

The two elements of the circuit are considered in detail in the later 
sections of this paper. Before this, their functioning is briefly described. 


(a) Amplification by the jet 

It is well known that, if one fluid moves over another with uniform 
velocity, the common surface may be unstable to certain disturbances 
(see, for example, Lamb (1932, p.373)). Particular examples are the 
generation of waves by wind and the flapping of a flag. If the fluids are 
of different densities, only disturbances of small wavelength (i.e. high 
frequency) are unstable, but, if the fluids have the same density, all 
disturbances are unstable. Surface tension stabilizes small wavelengths 
so that, with a suitable choice of parameters including the densities of the 
fluids, both large and small wavelengths may be stable, leaving only a 
narrow band of frequencies unstable. One frequency in this band will 
eventually predominate. Similar results apply to a circular jet. If it 
moves through a fluid of the same density, disturbances of all wavelengths 
are unstable as is shown by the breaking up of an open jet about one diameter 
downstream of the nozzle. A slotted wall may be regarded as a stabilizing 
influence, restraining the radial motion of the jet boundary. ‘The exact 
mechanism by which the slotted wall restrains the motion is not clear but 
it will be shown that both large and small wavelengths may be stable. The 
measured amplification is considerably less than that calculated for an 
open jet, thus indicating that the slotted wall has a considerable stabilizing 
influence even although complete stability is not achieved. 


(b) The return path 

At entry to the diffuser, the jet contains more oscillatory energy in a 
certain frequency band than when it left the nozzle. If a return path exists 
such that energy can be fed back to reinforce the original disturbance for 
one or more of the frequencies in this band, then large oscillations will occur 
at these frequencies. One such path, found in the wind tunnel, is formed 
by the reservoir. Ignoring turbulent mixing, there is no interchange of 
particles between the jet and the reservoir and only particles which have 
issued from the nozzle enter the diffuser. ‘The reservoir thus contains 
a constant mass of fluid whose volume fluctuations due to the oscillations 


of the jet boundary set up pressure fluctuations, 
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How these act on the jet leaving the nozzle is not clear; it is possible 
that they cause the nozzle to vibrate radially. ‘lo ensure a phase shift of 
an integral number of cycles round the loop, it is necessary that the jet 
fluctuations should produce the maximum pressure in the reservoir when 
the jet diameter is a minimum at the nozzle, giving the result that the length 
of the working section must be an integral number of wavelengths less one 
quarter, a result borne out by experiment. | It appears that this mechanism 
still operates in the water tunnel, but there is now another return path. 
The water tunnel differs from the wind tunnel in being of closed return 
type, and pressure fluctuations are transmitted round the return path of 
the circuit. ‘This results in certain frequencies, approximately independent 
of water speed, being amplified. 


3. JET INSTABILITY 


Derivation of the basic equations 

The cross-section of the tunnel working section is circular. The problem 
considered is the instability of a circular jet of incompressible fluid, of 
radius a, moving with velocity U through an infinite extent of the same fluid 
at rest. The disturbances whose stability is considered are axisymmetric, 
since only axisymmetric disturbances can operate the feed-back mechanism. 
(It is assumed that the effect of the slots and gaps of the slotted wall can be 
averaged round the circumference.) ‘The analysis has been carried out for 
similar situations by various authors (e.g. Rayleigh 1879), and so is given 
here only in outline. 

The two coordinates of a point are taken as 7, the radial distance from 
the axis, and g, the axial distance from the (arbitrary) origin. At any instant, 
the common surface of the jet and surrounding fluid is at 

r=ary, 
where 7 is a function of z and ¢, the time. It is assumed that 7 is of the 
first order of smallness; in particular, boundary conditions may be satisfied 
at r = a rather than at r = a+7. 

Inside the jet, the velocity potential takes the form 

® = Uz+4, 
and outside @ = ¢’, 
where ¢ and ¢’ are the small disturbance potentials associated with 7. 
These potentials satisfy Laplace’s equation. 

The observed disturbance is a progressive wave moving downstream. 

Corresponding to the surface disturbance 
= no exp[i(ot —kz)], 
we find the solutions 
d = Al, (kr)exp[i(ot —kz)], 
¢' = A’K,(kr)exp[i(ot — kz)], 


having the correct behaviour at r = 0 and r = ~« respectively, where J, 
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and K, are the modified Bessel functions of order zero in the notation of 
Watson (1944, pp. 77 & 78). (Although the reservoir is not of infinite 
diameter, calculations show that the simple expression for ¢’ gives numerical 
results of sufficient accuracy.) 

Of the three boundary conditions to be applied at r =a, two are 
obtained by relating the radial component of particle velocity to the motion 
of the jet boundary. ‘These give 


an _ ag’ 

ot Or ’ 
On U on _ of 
ot oz er’ 


the extra term in the second equation arising from the translational velocity 
of particles inside the jet. ‘he third condition relates the pressure at the 
boundary inside and outside the jet. By a mechanism which is discussed 
later, the slotted wall section is assumed to impose a pressure difference 
across the common surface, so that 


od U Ca og’ - 
ot Og ot p f 


where p is the density of the fluid, and P = P, exp[i(ot — kz)] is the pressure 
ditference, taken to be positive when the pressure is greater inside the jet. 
The boundary conditions give 
tony —tkU ny = ARI,(ka) = AkI,(ka), 
tony = A’RK/ (ka) A‘kK, (ka), 
ig Al,(ka) —tokRU AI, (ka) = to A’ K, (ka) — Py/p; 
and elimination of 4 and A’ from these gives 


os Oe 2g SO 
kI,(ka) kK,(ka) pn 
\ suitable non-dimensional form of this equation is obtained by the 
substitutions 
t-2., «See. é.29 2S 
kl ,(ka)K,(Ra) kal,(ka) U? py,’ 


which reduce it to 


I(ka) a P, 


7? 





(X —1)?+aX? = 8B. (1) 
The quantity z is a function of ka only, which increases from zero at ka = 0 
to unity as ka-> x. It is plotted in figure 3. 

‘There are two ways of investigating stability. Inthe situation considered, 
vis real and, if (1) gives complex values of YY, then & is complex and stability 
depends on the sign of the imaginary part. Although this represents what 
happens in practice, it is inconvenient to use this approach as x and f are 
complicated functions of the complex variable &. 

‘The usual way of investigating stability is to consider an initial sinusoidal 


surface disturbance and to investigate its behaviour at subsequent times. 
Then & is taken to be real and (1) is used to determine X as a function of k 
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for all values of k from zero to infinity. ‘The two cases are related by the 
assumption that the amplification per cycle in the latter case equals the 
amplification per wavelength in the former. Although this is not exactly 
true, it is approximately true provided that the amplification is not too large. 
Let a typical root of (1) be ¥ = X,+72X,, where X, and _Y, are real; 
then o = o, +10, = RU(X,+7X,). All the unknowns », ¢, ¢’ and P contain 
a real exponential factor exp(—o,f) and the motion is therefore unstable 
when o, < 0 or X, < 0. 
x 
10 








ala. ’ 4, 
re) ' 2 3 4 
Figure 3. Plot of « ws ka. 





The measurements taken at a given speed are the frequency f and 
wavelength A of the disturbance, and also ¢, the amplification in a given 
length /. It follows immediately that 

k= 2n/d, o, = 2nf, (2) 
and, on the assumption that the amplification per wavelength equals the 
amplification per cycle, 


o>) 
— 


e = exp(—Akla,/o,). (: 
The ratio of wave velocity to main jet velocity is 

fa[U =0,/kU = X,. 
Since x is a function of ka (i.e. 27a/A), x is known from the wavelength. 
By means of (2) and (3), k, o, and o, are known and hence X. ‘Thus the 
entire left-hand side of (1) may be determined frem the measurements, 
so that the effect of a particular slotted wall configuration can be assessed. 


The open jet 
For an open jet, there is no stabilizing mechanism and 8 = 0. ‘Then (1) 
becomes 
(X—1)?+ aX? = 0, 
F.M. fy 
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as was shown by Rayleigh (1879). Since x > 0, this admits only the complex 
solutions 
L+1a)? 





’ 


l+ 
one of which is unstable. The ratio of wave velocity to main jet velocity is 


: 1 
X eet ae 4 
1” 1+¢ 4) 
so that this ratio always lies between 0-5 and 1. The amplification in 
length / is 
e = exp(k/z!”). 


-_— 
wn 
~~ 


Particular forms of stabilization 

To discuss stabilization, it is necessary to consider the term f in (1). 
This is directly related to P, the pressure difference across the slotted wall, 
which we postulate to be the mechanism by which the flow is stabilized. 

The slotted wall may possibly set up a pressure difference in two ways. 
‘lo some extent, it may move in and out with the jet boundary. ‘Then it 
behaves as a skin surrounding the jet and its elastic properties give rise to 
the pressure discontinuity. Pursuing this analogy, one might expect a 
term due to surface tension; it would be expected that the bars bending 
and being pushed against their supports would also give rise to pressure 
discontinuities of this nature. 

The other function of the slotted wall is to permit flow through the 
gaps. ‘There will then be a pressure difference necessary to overcome the 
resistance to flow. ‘The total resulting pressure difference is a complicated 
non-linear function of 7 and its derivatives, involving, among other things, 
an unknown relationship between the amplitudes and phases of the motion 
of the jet boundary and the bars. For small disturbances, this function 
may be linearized and, although magnitudes cannot be determined, the form 
of several particular terms in the linear expression can be suggested. The 
stiffness of the bars on their supports gives a restoring force dependent on 
the displacement of the bars and so of the jet boundary. Bending of the 
bars gives a restoring force proportional to the fourth derivative in the 
axial direction of displacement. Resistance to fluid forced through the 
gaps gives a term dependent on the outward velocity of the jet boundary. 
It has been pointed out by Rayleigh (1879) that surface tension if it exists 
gives rise to two terms, a restoring force proportional to the curvature in 
the axial direction and a force proportional to displacement which tends 
to produce divergence from the equilibrium position. Thus we might 
expect P to contain terms in some or all of 7, 6?y/éz", é4y/éz* and én/ct. 

Substitution of the sinusoidal oscillations then gives P, in the form 


Pon = P, (ka) +1UP,(ka)X, 


where it is to be expected from the foregoing considerations that P, — const. 
as ka->(0), P;-> o at least as rapidly as (ka)? as ka— ox, and P, x ka, 
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Then 
I,(ka) a si ‘ 
B= —*—— — {P, (ka) +iUP,(ka)X} 
hal,(ka) pU? (ka) +71UP,(ka) 
= 2iyX +4, (6) 


say, where y and 6 are positive functions of ka. Hence (1) becomes 
X7(1+a)—2X(1+iy)+1-6 = 0, 
the roots of which may be expressed in the form 


(Z+1)(GY +1) 





X= ? / 
— (7) 
where 
Y2-Z? = y?+f, ) 
YZ =, Y,Z>0,} (8) 


From (7) we now have 
, 
l+o 
An unstable solution only exists if X, < 0. Thus the upper sign always 
gives a stable solution and the lower sign gives an unstable solution only 
if Z <1. 
For an unstable solution, the ratio of wave velocity to main jet velocity 
is X, = (1—Z)/(1+«), so that 
0< X, < 1/(1+2), (9) 
and the amplification in length / is 


e = exp(R/))). (10) 


‘lo investigate the range of instability, it is simplest to consider the 


. eu te 





(y*,¢)-plane. Z is constant along the straight line 

y*(Z-* -1)-C = 2’, 
and, as is seen in figure 4, the region 0 < Z < 1 becomes the region € > —1 
for all y?.. The special case y = 0 gives stability if ¢ < 0 and instability if 
€>0. The amplification depends on Y, which is constant along the 
straight line 

y(V2+1)+0= 9%, 
The line of neutral stability € = —1 corresponds to 6 = 1. 

Several general conclusions can be drawn from the form of £ given in (6), 
making use of properties of the modified Bessel functions given by Watson 
(1944). 

(1) For large wavelengths, i.e. Ra small, 

y~ Cha 4ol ’ 
where C is defined by P,(ka) = (C/a)ka, 
8 ~ P,(0)a/2pl?, 
and 4 ~ }(ka)* log(2/ka). 


Disturbances are stable if P,(0) > 2pU,a. 
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(ii) For small wavelengths, i.e. ka large, 
y ~ C/2pU, 
6 ~ Cha/pU?, 
where C” is defined by P,(ka) ~ (C’/a)(ka)?, and C’ > const. > 0 or C’ > 2 
as ka-> x, and x ~ 1. Disturbances are stable if C’ > pU?/ka, which is 
certainly so for large enough values of ka. 





Z = CONSTANT 


- Y¥ = CONSTANT 
4 


dl Madea Sia’ 























Figure 4. Region of instability in the (y*, {)-plane (¢ > — 1). 


(iii) Even if large and small wavelengths are stable, an intermediate 
band may not be. Figure 5 shows a particular case where P, (ka) takes the 
form 4 + B(ka)*. 

(iv) Increasing the speed of the jet U increases the range of instability, 
as is also illustrated in figure 5. It is seen that increasing the speed by 10°, 
in the example plotted makes all large wavelengths unstable as well as 
decreasing the lower limit of unstable wavelengths, and decreasing the 
speed by 10°,, stabilizes all wavelengths. 

(v) Damping increases the instability of the flow. For a constant 
wavelength, increasing y increases Y and so increases both the amplification 
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per unit length and the amplification per cycle. ‘Too much should not be 
made of this apparently paradoxical result. In the first place, it may be 
physically unrealistic to suppose that y can be altered without altering ¢. 
In the second place, the most relevant measure of instability in this theory 
is the rate of increase of the disturbance with time. It is seen that the 
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Figure 5. Theoretical example of the dependence of stability on wavelength and 
main jet velocity. 


disturbance increases with time as exp{kU Yt(1 — Z)/(1+«)}; this increases 
with y only if Y* < Z, which is only so inside the curve ( = y?*—y43—,? 
shown in figure 4. Nevertheless, it is clear that, in some circumstances, 
the statement is correct, though it should be noted that change of damping 
cannot destabilize a stable flow. 


4. ‘THE RETURN PATH 
Quite a broad band of frequencies may be amplified by the jet (as 
demonstrated in figure 5), and it follows that many paths may exist returning 
energy from the downstream end of the jet to the nozzle and providing the 
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correct phase relationship for one or more frequency in the band. However, 
in most cases, the efficiency of the return path is likely to be so low as to 
reduce the gain round the loop to less than unity, thus eliminating any 
instability in this loop. It is not possible to consider return paths in any 
generality, and any new tunnel will present its own problems. ‘Two particular 
return paths have been identified in the A.R.L. wind and water tunnels, 
and their occurrence and ettect on the tunnels will be described. 


Role of the reservou 

goth return paths have one element in common. ‘The boundary of the 
jet consists of oscillating particle paths leaving the nozzle and re-attaching 
at some point near entry to the diffuser. It thus encloses a constant mass of 
fluid in the reservoir which, in the theoretical inviscid, non-oscillating 
situation would be at rest, separated from the jet by a free streamline. 
(The effect of viscosity, together with turbulent mixing, is to thicken the 
jet boundary and to impart momentum in the direction of the jet to the 
neighbouring fluid in the reservoir. By continuity, this fluid is forced to 
return along the outer wall of the reservoir so that the flow in the reservoir 
takes the form of a general swirl set up by an annular ring of vorticity.) 
The oscillating jet boundary produces a periodic volume fluctuation of the 
fluid in the reservoir and so an oscillatory pressure. ‘This explanation 
introduces the compressibility of the fluid. However, this does not have 
much effect on the amplitude and phase distribution of pressure in the 
reservoir, as is easily verified by considering the reservoir as a tube closed 
at the nozzle end and energized by the vibrations of the downstream end 
of the jet boundary. In all practical cases considered, the tube 1s well less 
than a quarter wavelength long, and so the oscillating pressure does not 
vary greatly in amplitude or phase throughout. ‘This pressure fluctuation 
must also exist throughout the jet, giving rise to a fluctuation in density ; 
energy must be propagated out of the working section both up and 
downstream. Rough calculations indicate that, in all cases considered, 
these pressure fluctuations should completely mask the fluctuations 
in the jet which are calculable from the theory of the preceding section 
in the form of amplified progressive waves. 


Return path through the reservoir 

One return path utilizes the potential energy of these pressure fluctuations 
in the reservoir. How energy is fed back to the jet is not clear, but it seems 
probable that the nozzle is set into vibration by the fluctuations. On this 
assumption it is possible to calculate the preferred frequencies of this 
circuit. 


The displacement of the surface of the jet may be expressed as 


n = ny exp[i(ot —kz)], 


where k is complex. The change in volume of the reservoir is V = I, e' 
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with 
el 
V. = 2man, | e-'** dz 
/0 


_ 27AN) poke: 7 pa 
= ——— {e'sink,1+71(e*:'cosk,l—1)}, 

k, +tk, 
where k = k, + thy, ky, ky > 0. 

When V takes its maximum value, the reservoir volume is least and so 
the pressure is greatest. ‘This pressure acts on either the jet or nozzle at 
the upstream end and it is assumed that at this point the jet diameter is 
then a minimum. ‘Thus V and 7 must be an amount z out of phase, so that 
Vy/y) is real and negative. ‘This condition is only satisfied if 

e*'sink,! — ek*!cosk,1—1 








ke, k <& 
Experimentally, it is known that e“’ 5 1 and k, > k,. Thus sink,/ = —1, 
and 
k,l = 2nr+ 30/2 + &, (11) 
where 7 is an integer. 
Since € is small, 
es! = Egksl_ 
~~ = 


so that 

£ = e~*e!— fh, /k,. 
Experimental results indicate that e~":’ and k,/k, are comparable in 
magnitude, so that € = 0. Thus, from (11), 


lid = n+ 3. (12) 


Return path round the tunnel return circuit 

The other return path utilizes the energy propagated out of the working 
section. It is applicable only to a closed return tunnel when the energy 
may move round the circuit with relatively little loss. ‘The preferred 
frequencies are then those for which the time of transit round the tunnel 
circuit is an integral number of periods. 


5. EXPERIMENTAL WORK 

Description of tunnels 

(a) Water tunnels 

The 30 in. tunnel is of closed return type as shown in figure 1, and, 
apart from the working section, is only noteworthy here for its ‘resorber ’. 
This feature, whose use in water tunnels originated in the U.S.A., is a 
portion of the tunnel circuit so located that water is subjected to pressure 
for long enough to redissolve any air bubbles which may have formed during 
cavitating flow in the working section. In the present work, its relevance 
is that it provides an air trap in the circuit halfway round the resorber which 
materially affects the vibrations observed. 
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‘The working section is six diameters long, and the reservoir diameter 
is two and a half times the jet diameter. ‘The slotted wall consists of sixteen 
rectangular bars, the gaps between bars comprising 20°, of the total 
boundary. Further details of the tunnel are given by Lever et al. (1957). 
‘he tunnel was equipped for these tests with static pressure holes in the 
wall of the reservoir, diffuser and contraction, and Pitot-static tubes in the 
main jet. 

A small amount of work was also done in the A.R.L. 12 in. tunnel. 
‘This was originally used as a prototype for the 30 in. tunnel, and is a 
scaled-down version of it in all important respects, except that it has no 


res¢ rl eT. 


(6) Wind tunnel 

The 7 in. diameter wind tunnel was designed specitically for the purpose 
of investigating slotted wall working sections. It is a blow-down tunnel, 
ending in an open jet to which can be attached any appropriate slotted wall 
unit. For these investigations, the dimensions of the water tunnel working 
section and diffuser were scaled down to fit on to the contraction and, after 
a reasonable length of diffuser, terminated abruptly, giving an open return 
circuit. ‘The slotted wall was made of Perspex, the outer reservoir wall 
and the ditfuser of acetate sheet, and the nozzle was a Dural extension to 
the wooden contraction. Some modifications were introduced during 
tests; in particular, the original contraction from the reservoir to the diffuser 
was replaced by a ring (see figure 10). ‘hese were shaped in wax with an 
appropriate stiff backing. Measurements were made with hot-wire 
anemometers, measuring the spectrum of the longitudinal component of 
turbulence with a frequency analyser whose bandwidth was about 1", of 
the mid-band frequency. ‘The sensitivity of the hot wires was not known 
with accuracy, and little use has been made of absolute values. Pressure 
measurements were also made, a static tube being joined to a microphone. 
Again absolute values were not known. 


Fet instability 

The first experiments in the small wind tunnel made use of smoke, 
and demonstrated that the disturbances were axisymmetric. Using two 
hot wires, it was possible to demonstrate by comparing phases that the 
disturbances were progressive waves moving downstream. ‘lo determine 
wavelengths, it was found simpler to add two hot-wire outputs which tend 
to cancel when the hot wires are spaced a half wavelength apart. Details 
of these results are presented later; we note here that these measurements 
indicate the ratio of wave velocity to main jet velocity to be approximately 
(1+)! so that, by (9), 7 < 1. 

Removing the outer reservoir wall in this tunnel eliminates the instability, 
and the spectrum of the flow issuing from the nozzle is fairly ‘white’. 


Comparison with spectra measured downstream gives the amplification 








aes 
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along the working section as a function of frequency. A typical example, 
at a tunnel speed of 75 ft./sec, is given in figure 6 and indicates that at this 
speed frequencies between 20 and 60 c/s are amplified with a maximum 
amplification along the length of working section of 6:6 at 40 c/s. If the 
wavelength is calculated on the assumption that Z < 1, it is found, using (10), 
that Y < x'* < Lalso. Thus, from (8), y <1, ¢| « 1, andso6 = «/(1+z2). 
It is of interest to compare the results with those calculated for an open 
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Figure 6. Measured amplification along the jet as a function of frequency for the 
7 in. wind tunnel at 75 ft./sec. 


jet from (4) and (5). The amplification along the working section for an 


open jet is found to be a rapidly increasing function of frequency, which at 
75 ft./sec, takes the values 13-5 at 20 c/s and 360 at 30 c/s. The effectiveness 
of the slotted wall in stabilizing the flow is amply demonstrated. The very 
much less extensive measurements for the 30 in. water tunnel also confirm 
the tendency towards stabilization. 
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Return path 

(a) Wind tunnel 

The wind tunnel is of open return type. ‘Thus disturbances can only 
return through the reservoir and the preferred disturbances are those with 
an integral number of wavelengths less one quarter in the length of the 
working section. 

Since frequency is more easily measured than wavelength, figure 7 (a) 
shows the observed peak frequency or frequencies as a function of tunnel 
speed. ‘To calculate the expected frequencies, it is assumed that Z < 1 
so that, from (9), the ratio of wave velocity to main jet velocity is (1+)~’. 
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Figure 7. Dependence of instability on tunnel speed for the 7 in. wind tunnel: 
(a) preferred frequencies and wavelengths; (b) turbulence level at entry to 
diffuser. 

For each possible wavelength, frequency is then proportional to tunnel 

speed, and the corresponding straight lines are shown in figure 7(a). ‘The 

theoretical wavelengths and some of the measured wavelengths are also 
shown. This leaves no room for doubt that the basic mechanism outlined 


above is correct and, as previously stated, that the assumption Z<1is 
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justified. ‘lhe most interesting feature is the jump in frequency at 75 ft./sec, 
corresponding to a change from two and three-quarters to one and three- 
quarters wavelengths. Figure 7 (6) shows the turbulence level at the diffuser 
end of the working section as a function of speed. Frequencies in the range 
about 30 to 40c/s only are accepted, presumably due to a mechanical 
resonance associated with this particular tunnel. Below 75 ft./sec, the 
mode with m = 1 (see (12)) is below this frequency range and, above this 
speed, the mode with nm = 2 is above it. At 75 ft./sec, both modes are 
amplified simultaneously and figure 8 (a) shows the spectrum at this speed 
together with that obtained at the same speed with the outer reservoir 
wall removed. In figure 8 (4), a short length of trace is shown taken from 
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Figure 8. Typical records of turbulence at entry to diffuser at 75 ft./sec. in the 7 in. 
wind tunnel: (a) spectra (i) with outer reservoir wall in place, (11) with outer 
reservoir wall removed ; (4) cathode-ray oscilloscope record of hot-wire 
output with outer reservoir wall in place. 


a cathode-ray oscilloscope record of the hot-wire output showing the beats 
produced by the addition of the two frequencies. Since jet instability 
increases with speed, it is surmised that modes with n = 3 and above are 
not observed since the gain along the jet is too small at the relevant low 
speeds. The mode with n=0 is never observed since the relevant 
frequencies only arise above the top speed of the tunnel. 
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(6) Water tunnels 

The problems of measuring fluctuating pressures and velocities in 
water are very much greater than those in air. Consequently, although 
this investigation was aimed at improving the 30 in. water tunnel, the 
phenomena have been less comprehensively investigated there. ‘The 121n. 
water tunnel fares even worse, as its smaller dimensions and the higher 
frequencies arising increase the difficulties. All that can be said of the 
12 in. tunnel is that, although the oscillations were not initially noticed, 
they were observed when looked for later, and the frequencies observed 
confirm the existence of the feed-back loop of the wind tunnel. 

For the 30 in. tunnel, continuous vibrations such as are obtained in the 
wind tunnel were rarely obtained, spasmodic bursts being more usually 
encountered. Frequencies were measured by counting cycles on the 
oscilloscope traces, and in some cases beats were observed. Although the 
beat frequency could be estimated, an ambiguity arises as to whether the 
second frequency is above or below the parent frequency. In the presen- 
tation given here, this ambiguity has been resolved simply to fit in best 
with the general picture presented. 

It was discovered—inadvertently in the first instance—that the results 
were dependent on the amount of air trapped in the resorber. With no 
air in the resorber, extremely large disturbances were obtained at tunnel 
speeds between 45 and 60 ft. sec, the shake of the tunnel both being visible 
and imparting considerable vibrations to the floor. (In the working section, 
pressure fluctuations up to 9 |b. 1n.2 were observed.) ‘The presence of 
an air bubble reduced the disturbances to levels which were subjectively 
acceptable at all speeds. It was therefore surmised that, with no air bubble 
in the resorber, the return path was by way of the closed return circuit, 
and the air bubble, when present, acted as a pressure release surface to 
eliminate this return path. 

‘This surmise is not supported by the frequency measurements, which 
are presented in figure 9. With air in the resorber (figure 9(a)), the 
observed frequency is almost independent of water speed being about 
6-5 c/s with most observations lying between about 6-2 and 6:8 c/s. With 
the resorber vented (figure 9(6), which also includes some anomalous 
results with air in the resorber), the observed frequencies seem to fall 
into several groups, each increasing with speed although not so rapidly as 
predicted by theory. Again considerable scatter is observed; this ts 
known to be genuine as different bursts in the same run can reveal marked 
variations both in frequency and in the strength of observed beats. ‘lhe 
observed scatter may be due to the practical difficulty of controlling the 
amount of air in the resorber. 

A theoretical estimate of the time for an acoustic pulse to traverse the 
return circuit of the tunnel gives 153-2 milliseconds, so that a preferred 
frequency of 6°53.¢s would be expected for this path. It is clear that, 
with air in the resorber, this path dominates. With the resorber vented, 


the amplitude increases and this appears to stimulate the alternative return 
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path through the reservoir. Then, above 40 ft./sec, the two return paths 
appear to operate in parallel, frequencies between those predicted for the 
return paths separately being obtained. At the top of the speed range, 
one or two results suggest that the mode for which m = 3 may dominate and 
impose its frequency on the entire motion. 
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Figure 9. Dependence of preferred frequencies on tunnel speed for the 30 in. water 
tunnel: (a) with air in resorber; (4) with resorber vented. 


This whole picture suggests that complicated non-linear effects are 
present (for example, a threshold of instability and frequency entrainment) 
which are not taken into account by the theory. Nevertheless, it is considered 
that the evidence does confirm the existence of the two return paths, 
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In both the 30 in. water tunnel and the wind tunnel, the conclusion 
that the pressure is constant in amplitude and in phase throughout the 
reservoir and working section is borne out by measurement. 


6. METHODS OF ELIMINATING INSTABILITY 

The purpose of this paper is to give an explanation of the phenomena 
observed in slotted wall tunnels and to illustrate this with the evidence 
gathered to date. However, since the understanding of the cause of the 
oscillations indicates several lines of approach to eliminate them, it is of 
interest to list these briefly. The degree of success achieved with the 
various suggested schemes is critically dependent on unknown parameters 
of the particular tunnel, especially in the dependence of amplitude on the 
degree of non-linearity of the circuit as well as the gain. For this reason, 
numerical results are only given to illustrate the discussion, and no great 
detail is entered into. The particular measures eventually adopted to 
cure the water tunnel will be reported elsewhere when the investigation 
is completed. 

The first place in which one looks for improvement is the jet itself, the 
amplifier element which is essential for any instability. ‘The stiffness of the 
slotted wall can be increased either by stiffening all members or by reducing 
the gap width. In the wind tunnel, reducing the total gap width from 20", 
of the circumference to 7°,, has produced a worthwhile reduction in 
amplitude. ‘The optimum gap width is likely to be determined by the failure 
of the section to operate as a slotted wall section, giving model corrections 
comparable to those of a closed jet. Another remedy lies in the fact that, since 
the length of the working section is of the order of only two or three wave- 
lengths, end effects may be expected to be significant. ‘This is exemplified 
in figure 10 which illustrates three different types of entry to the diffuser 
tested in the water tunnel: viz. a contraction, a ring and a splitter ring. 
It has been found that the splitter ring is very successful while the ring gives 
results rather worse than the contraction. ‘The reason for this is not clear. It 
is conjectured that the jet boundary at re-attachment to the wall is effectively 
moored by the splitter ring whereas the contraction gives too much freedom 
of movement to the re-attachment point and the ring to the direction of 
re-attachment. Whatever the reason, this confirms the point that end 
effects may be utilized to reduce instability. 

Common to both return paths so far observed is the production of large 
pressure fluctuations in the reservoir. ‘The absolute pressure fluctuations 
are proportional to the proportional changes in volume of the reservoir, 
and thus may be reduced by increasing the size of the reservoir. ‘This has 
not been tested but, as shown in figure 8 (a), removing the reservoir wall 
(thus making its size infinite) completely eliminates the oscillations. 
A pressure release surface in the reservoir would also have a similar effect. 


This involves considerable constructional difficulties in the water tunnel; 
in the wind tunnel the effect has been simulated by drilling holes in the 
outer reservoir wall, Removal of about }°, of the total wall area was sufficient 











Instability in slotted wall tunnels 303 


to reduce the fluctuations to one-tenth of their value. For the two particular 
return paths, the mechanism by which the pressure fluctuations in the 
reservoir feed energy back to the nozzle is not yet understood. Accordingly, 
no remedies can be offered. The efficiency of the return path through the 
remainder of the closed return circuit can be reduced by providing a pressure 
release element, as has already been shown. 


(a) Le. aS 




















Figure 10. Types of entry to the diffuser tested in the 30 in. water tunnel: 
(a) contraction; (4) ring; (c) splitter ring. 


This section thus provides further evidence that the diagnosis of this 
paper is correct by showing that deductions from the theory are in accord 
with experiment. 


The authors are indebted to many of their colleagues and particularly 
to Dr H. Ritter, who was responsible for the tests in the water tunnel and 
also for several valuable comments. 


APPENDIX. "THE EFFECT OF COMPRESSIBILITY ON JET INSTABILITY 
In view of the use of the slotted wall section for wind tunnels under 
transonic conditions, it is of interest to consider the effect of compressibility 
on jet instability. 
As in the incompressible case, 


cM 0 
; w= —, where ® = Uz+d¢. 


b= 





~ 


or (ores 
For small disturbances, we put p = py+p’, p = pyt+c*p’, where c is the 
velocity of sound, and we neglect second-order terms in ¢ and p’. The 








304 j. L. King, P. Boyle and 7. B. Ogle 


equations of continuity and motion then give 
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The substituting of a solution of the form 


$ = $o(7 )exp[i(ot — kz)] 
gives an equation for ¢y: 
d°hy | 1 dbo _ ,, , : 
——~ + -— —k{1—M*(X—1)}*}4, = 0, 
dr? r dr ( Yio 
where M is the free stream Mach number U/c and X = a/kU as before. 
Using the boundary conditions, we can derive an equation for X in the same 
form as before: 
(X—1)?+a’'X2 = fp’, (1’) 
where 
sil I,(k'a)K,(k’a) 
I)(k'a)K,(h'a) 
_ Ral,(k’'a) 1 a@ Py 
1)(k‘a) kia? U2 Po To ; 
2= k*{1]- M*(X-1)*}. 
The appearance of X in x’ and ’ through its appearance in k’ immensely 
complicates (1) which is otherwise identical with (1). In particular, we 
are looking for complex values of X with the result that k’ is complex. 
When is small, we can proceed formally and expand in terms of WV. 
‘Taking the first correction term only, we put 
I,(k’a) = I,(ka) — §M?(X —1)?kal/ (ka) 
with similar expressions for Ky, /, and K,, so that 
x’ =%+M2(X-1)%, p’ =B+M{X-1)%8, 
where « and § are the values for the incompressible case, and 


(L(ka)  K,(ka)  I)(ka) K,(ka)) 








E ka. + — = Sees >y 
x (Ij(ka) K,(ka) I,(ka)  K,(ka); 
B shal Cha) I,(ka)| 

B \I,(ka) —‘1,(Ra)§’ 


in which well-known recurrence relations for Bessel functions have been 
used. If we replace XY by ¥ + .W?(X—1)?2X, where XY is now the incom- 
pressible value given by (7), it is found that, to first order, 
2NIN(1+%)—(14iy)} = —2X24 29X48 
= —@X?+(B/B)X2(1+2%)—-2X +1}, 


since y/y = 5/5 = B/B. 
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It has been seen in the main part of this paper that, in practice, Z < 1 
andy < Y<a!*<1. Then X = (1—7Y)/(1+z<) for an unstable solution, 
and 


-2iX(Y¥ +y) = 





ot — - 
— —_——. (1—2:Y— Y*)+ — 

(1 " a)? ( : ) A 
or, approximately, 


Re i = pee 1 ~(1+a) Bl. 
(l+a)? 2Y(1+a)? la B} 
The real part of X is positive, and the much greater imaginary part is 
negative. If it can be assumed that Y* < x’, which is reasonably certain 
in practice, then it is seen that the effect of compressibility is to increase 
the amplification of an unstable flow. 

The presence of a body in the tunnel affects the flow and the local Mach 
number may difier greatly from the free stream Mach number. Thus, 
although at low Mach numbers the presence of a body should only have 
a small effect on the stability of the tunnel, at higher Mach numbers the 
effect may be considerable, and a more elaborate theoretical attack does 
not seem justified. 
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SUMMARY 

In this paper the propagation of long waves (tides) into a canal 
is studied. ‘The canal is assumed to be rotating at constant angular 
velocity and the depth of the fluid is uniform. 

The rate of rotation can have a considerable effect on the 
amplitude of the unattenuated modes in the channel. ‘This is due 
partly to the modification of the known solution when the rotation 
is zero, and partly to the fact that even if there is only a single 
semi-infinite barrier unattenuated waves of a special type (Kelvin 
waves) may be propagated in the rotating system into what is 
normally called the ‘shadow’ region behind the barrier (Crease 
1956). 

The purpose of this theoretical investigation is to seek a partial 
explanation of the behaviour of tides and storm surges in the 
North Sea. For this reason two models of this region are discussed. 
In Example I the model is of two parallel semi-infinite barriers in 
the path of a plane progressive wave, and in Example II there are 
two parallel barriers, one semi-infinite and the other infinite. 

The ratio of the observed height of the semi-diurnal M, tide 
in the North Sea to the height of the incident tide from the 
Atlantic lies between the results predicted by the two models and 
is in closer agreement with the model of Example I. 


1. INTRODUCTION 

It has been shown in a previous paper (Crease 1956) that long waves 
incident normally on a semi-infinite barrier in a rotating system give rise 
to waves of Kelvin type which are propagated without attenuation into 
the region behind the barrier. ‘The direction of propagation is parallel 
to the barrier and the crest heights decrease exponentially in the direction 
normal to the barrier. When the ratio of the wave frequency o to the 
angular velocity w of the rotating system ts such that 2w o > 3/5, the amplitude 
of the Kelvin wave at the barrier is greater than that of the incident waves. 

This effect appears to explain qualitatively the formation of the semi- 
diurnal tides in the North Sea; for, to a first approximation, the British 
Isles may be regarded as a semi-infinite barrier across the path of the M, 


tide propagated in from the Atlantic, 
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In this paper the effect of the continental coast on the formation of the 
Kelvin wave is examined. ‘Two extreme cases are considered as an 
approximation to the continental boundary. In Example I (§3 to §5) the 
boundary is taken to be a second semi-infinite barrier parallel to the first 
and a distance 2a behind it. In Example II ($ 6 to § 11) it is assumed to 
be an infinite barrier distant 2a behind the semi-infinite barrier. 

In neither of these systems is there a barrier corresponding to the 
southern end of the North Sea. Taylor (1921) has considered the reflection 
of Kelvin waves at the end of a long gulf, but here we are only concerned 
with the transmission of waves into a channel from outside. ‘The possibility 
of resonance effects is thus ruled out. 

The problem when w = 0 is a familiar one in acoustics and scalar electro- 
magnetic theory, and Heins (1948) has solved the problem of an incident 
electric field parallel to the edges of the barriers. It is his method of 
solution, based on the Wiener-Hopf method, that is used in this paper. 
Vajnshtejn (1948) has also solved this problem and the acoustic problem, 
and investigates his solutions in more detail. 

In the case of a single barrier, rotation gives rise to the Kelvin wave 
mentioned above in addition to the usual diffracted wave diverging from 
the edge of the barrier. When there are two semi-infinite barriers the 
parts of the diffracted waves propagating into the channel may be regarded 
loosely as giving rise to the ‘channel’ modes of acoustics (w = 0) which, 
by the rotation of the system, become waves of Kelvin and Poincaré type 
(Proudman 1953, p. 265). In addition, the Kelvin waves arising from the 
barriers individually also propagate down the channel. ‘The amplitude of 
the waves in the channel will depend on its width, which in this paper is 
restricted to be less than half a wavelength. ‘This condition is satisfied for 
the North Sea. The observational data for the North Sea lie between 
the two cases, and is in better agreement with the model proposed as 
Example I. 


2. 'THE DIFFERENTIAL EQUATION FOR THE WAVE HEIGHT 
The equations governing the propagation of long waves in a system 
rotating at angular velocity w are, subject to certain approximations (cf. 
Lamb 1932, p. 318), 


ou fi oO 

eo oe 

ot ga? 

ov of 

= u = — g — 1 
ct f ° ay’ ( ) 
ou Ot 1 of } 

ee ae ee 

ox oy hot’ 


where u, v are the horizontal components of particle velocity parallel to x, y, 
h is the depth of water, and f = 2w. From equations (1) the equation for ¢ 
U2 
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is found to be 
ae See ee 
( “S ( a 5 1 ( a — 
ee ( a) (2) 
The boundary condition at a barrier is that the normal component of 


velocity is zero. When the barrier is along the x-axis, we have v = 0 and, 
in terms of ¢, 


— 
Ww 
— 


on the barrier. 
If the motion is periodic with time, we may replace ¢(x,¥,t) by 
C(x, y)exp(—tot), say. Equations (2) and (3) then become 


‘ar ay 
— + —- +e = 0, (4) 
ae oy 
where k? = (0? —f?)/gh, and 
oo if al . 
Saas (5) 
oy COX 


onthe barrier. It will be convenient to let 7f/o = tanzB, and to assume that 
k has a small positive imaginary part which may be taken to be zero when 
the solution of the problem has been obtained. 


EXAMPLE I. TWO SEMI-INFINITE BARRIERS 


3. THE INTEGRAL EQUATIONS 
Let the barriers be at y = +a, x > 0 (figure 1). By Green’s theorem 
(x, y) may be expressed as a contour integral round a contour I’ surrounding 
the point (x,y). ‘Thus, 


C(x, y) F | Gls Xo, Vo) = {( x0, Yo) — (Xr ¥o) — Gly; Xn do) | dso, 


0 Ong 

where 0/én, is differentiation along the outward normal to I’, and ds, is 
an element of I, which is taken to be a circle indented round the barriers 
and whose radius tends to infinity (figure 1). 

G is the free-space Green’s function given by (Morse & Feshbach 1953, 
p. 811) 

G(x, V3 Xo Vo) = FHA [R{(x — x9)? + (y —v9)?}2?], (6) 

where H is the usual notation for a Bessel function of the third kind. 

By using the boundary condition (5) and integrating by parts (Crease 
1956), it is found that 


P 1 Ge a se ae 
C(x,y) = — | < [f]{ cos78 — +siniB — ]G + 
i cos1B ! 9 CVo OX Yy=—a 


- { [I( cos — +sinizB a} | +¢, (7) 


q CVo OXy Ju,=a 





where [¢] is the limit as 6 +0 of C(x9, yp +8) — C(x, ¥p — 8), and it is assumed 
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that [C] is bounded at infinity. Here ¢; is the prescribed incident wave and 
is taken to be a plane wave at an angle 4, to the barrier. ‘Thus, 


¢; = exp7k(x cos 6) + vy sin 4). 
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Figure 1. 


The boundary condition is that 
cial ee 
cosiB — —siniB — }¢ 
OV Ox 


shall vanish on y = +a, x > 0. Thus, by operating on (7) with 
oe cao 
cosiB{ cosi8 — —sinip — 
oY Ce 


and setting y = +a, x > 0, we obtain the following pair of integral equations 
of Wiener-Hopf type for [¢(x,a)] and [¢(«, —a@)]. When these are solved, 


(7) gives C(x,v). Thus, on y = Fa, 


pce ee: 
219(%) = G,2(x) +4 (cosis — —sinifp ) . 
ae) 2 


Cy CG 


o rs A A 
{ . P ( a ” Oo ? 
| +t Iil% (cosie — +siniB ~)Gtwy: Xo» Vo) } 
iP cece . OV) OXp No a 


} | falsa)(co iB — +siniB 5 JOM V5 Bes Ve) + dx,, (8) 
OVo OXo Yo =a) 
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where 
f(x) = [C(x, ¥4)] v>0, 
= () Uy 
gi o(%) = tkeosif sin(#, — iP exp tk(x cos 4) | asin 4) x > 0, (9) 
= () =< U. 
gya(x) = 0 x >0, | 
and g,.(«) are defined by equations (12) for x < 0. These and following 


pairs of equations have been written in a condensed form. ‘The upper 
and lower signs in the equations refer to the first and second sufhces 
respectively (e.g. f, o(%) = [C(x, | a)] implies the two equations 
fi(x) = [¢(x, —a)], fo(x) = [€(x, @)]). 
Equations (8) may be added and subtracted; and since 
G(x, —a; X%),—@) = G(x, a; Xo, a) 
and G(x, —a; X»,a) = G(x, a; X%, —a) = G(x, 0; x», — 2a), 


we have 


~ 


ning sate Rn 
234(%) = g3.4(x)+( cos78B — —sinifp — } x 


OV Cx 
j ss hal , : 
| fsa(x){ cosiB — +siniB — }(G,, . »,+Gy, ~0) 4% 
J 2 oYVo OXy : ; 

(10) 
on y = 0, where g3, = g,+¥ etc. We may now take Fourier transforms 
of these equations, denoting transforms by corresponding capital letters, 
and obtain 

' a ee ae Sats 
G,,4(%) = O5 ,(%) + Fs ,(«){ cosi8 — —1xsin7B }{ cos78 — —1asinip } x 
oy CVo 
x {K (a, ¥,Vo)y =y,=0> K(a,¥, Yo)y =0,Y) = om * (11) 
where A(x, v,¥») is the transform of G(x,y;0,y)) and is given (Heins 
1948) by 


K(x, ¥, ¥o) = Saexp[?! y — v9) (Rk? — x”)! ?]/(k? — a?)P?. (12) 
Equations (11) may be conveniently written as 
G3.4(%) = O3.4(%) — F.4(%)A5.4(%), (13) 


in which A, ,() are explicitly 
K3,4(«%) = (x® — k® cos*78)(k? — x”)! ? exp ia(k? — x”)! i. in Jae — 4)1 2, 

(14) 
Also, from (9), 


O, ,(%) = 2k(x—R cos 6)! cos 78 sin(6, ia) sap (ak sinf,). (15) 


x —72S1N 


These equations are determinate, as all the transforms have a common 


strip of regularity $ Sk < Jia} < 0 (or < Sikcos6,! if this is 
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negative). The argument is the same as for the case of a single barrier 
(Crease 1956). The transforms K, ,(~) may be decomposed into the form 
K34(a) = Ky,()/Ky4(%), where A,,(«) are regular and free of zeros for 
JI \x} < F{k}, and Kj,(«) are regular and free of zeros for /(a} > — F(R}. 
Equations (13) may now be written as 
G34(a)A, (%) Os 4(%) | K, (a) - K, (A cos 65); 
Os (a) AS, (R cos 6) — F, (a) Ky (a). (16) 
In this form the leit-hand side of these equations are regular for 4 {x} >-4{R}, 
and the right-hand side for .7/x! < O(or < -*Aikcos6,' < 0). Thus each side 
of both equations (16) is an analytical continuation in their (overlapping) 
half-planes of entire functions F(x) respectively. Now A, ,(«) are 
O(x-!?) and K,,(«) are O(a!) as « — x (see below), and it follows 
from the known behaviour of the other functions that the left-hand side 
is O(z~1?) and the right-hand side is O(z!?). By an extension to Liouville’s 
theorem (‘Titchmarsh 1939, p. 85), the entire functions must be constants ; 
and, from the behaviour of the left-hand side as |x — ~, the constants are 
zero. Therefore it follows from equations (16) that 
F, ,(x) = O4.,(%) Kz, (R cos 0))/K, , (x). (17) 
It remains to determine the factorization of equations (14) explicitly. 


From (14), 





; Ks («) (x? — k? cos?78) . i.” 12 
K. x <=  Gemeceeueneaes = oe 7 —— k? = 1/2 5 
3(%) K, (x) (k2 — x2)!2 iexp| - ( e) 
f an-f — +" % tan (= ‘ 
; k— A k+ ] 
4a*(k? 77) 
| '- Gaye | 
"Thus 


4—kcosip 21. e+ a\12 
Kr (a) =14 ‘oar exp| as (hk? — x”)! 2 tan (7) vol) | x 


a 4a*k L’2 2iax 2iax 
TI ( (2n— a) (2n - 7 | so] (2n—1 = | 
4+kcosif 21 k— o\12 7 
— a+ Roos) exp| (k? — «”)!2 tan eat Xo(z | x 


K; (2) (k+2)!? 


4k? Le 2iax 21ax, 
aS SS ex EES - 19 
IT (1 (2n ae (2n—1 = | rl Ge 7 | om) 


where x (%) is an entire function introduced to make A, («) and K, (z) 
algebraic as x —- x. Asymptotically we find that as x — x, Fia} < F{R} 
(cf. Heins 1948) 


K; (a) ~ — (a/2)"* exp ak y +log = 


bo} 
a | 
w= | > 
= 
o< 
— 
R 


where y is Euler’s constant. 
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‘Theretore, 1f 





1a4. im 
Ky (%) — | 1—y+log (20) 
; TT | 2ak : 
then Ky, (a) ~ —(a/2)!* as jaj >a, F{a} < S{k}; and similarly 
Ks (a) ~ (2/x)"2 as ja\ > ow, Ff{a} > Zik}. Similarly, for A,(z), 


oi Ve 12 
K, (x) a(z—kRcos iB)exp| (k? ad 2 tan (G—) + (2) | x 
7 — % 
a*k*\12 1aqa7 1a ‘ 
i> exp( — —}, (21) 
: n> 7? nt nt 


| 21 _ y\ 1/2 
K; (x) (x+kcos esp = (Rk z2)\2 tan (G=) xl) | ; 


with 


t 
| 


The asymptotic behaviour of A, (a) and Ay (x) is then 
Ky (x) ~ Sa(iax/7)!?, KZ («) ~ (—a/iax)!?. 

‘The solution of the integral equations is now given by the Fourier 
transforms of equations (17), but these transforms will not be needed 
explicitly as € may be expressed directly in terms of F'3 4(). 

4. "THE WAVE HEIGHT € 
‘The integrals in (6) are convolutions and may be rewritten as 


; I ial ae , eg ee 
C(x, y)= earner | | Fite) (cosis ~ 2siniB) K(x, ¥,99)} 
477 COStp . ie \ 


\ Oo ja~ =< 


if e — os 
+ F,(%)< (cos iB — —ta sini) K(ay.y0)} Jess dx + 








OVo J 

texpthk(xcos4,+y sin 4) (24) 
»} < 0). Particular interest lies in the 
wave amplitude at large distances from the edges of the barriers. This 
is given by the residues at the poles of the integrand; the branch line 
integrals obtained by deforming the contour only give a contribution which 
tends to zero at large distances from the edges. The poles are at « = kos, 
kcosif, and the zeros of the infinite products of (18) and (21). ‘These 
latter zeros give rise to terms decaying exponentially to zero in the 
y-direction, provided 0 < 2ak < 7. 


where —-4{k! < « < 0 (or <.4Aikeos6 


j 


We consider first the pole at x = kcos@. If x < 0, the contour of 
equation (24) may be closed below, and there is no contribution of the 
pole to the integral. If x > 0, the residue is found to be 

expik(xcos@,+vsin@) for vy > —a, 
and sin(@, —18) 





exp ik{x cos 6) — (vy + 2a)sin 45} (25) 


sin(8, +78) 
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for y < —a. Thus, these terms together with the incident wave represent 
the reflection pattern obtained by a geometrical construction (apart from 
a phase change in the reflected wave due to rotational effects). 
Next we consider the residue at the pole « = kcosif. There is again 
no contribution for x <0. For x > 0, y >a, the residue is found to be 
4cosiB - 
sin(6, +78) 


1ak ies it \) 
x exp| F . 4, sin 6, +78 siniB + (cosiB — cos a)(1 —y+log . a x 
{ 2ak J 


x [cos }4) sin 378 cos'*(ak sin 9, )cos! *(ak sin iB )exp ty, + 





+ S1{sin 6) sin78 sin(@k sin 8))sin(ak sin 78 )}!? exp ix, | x 
x expikixcosiB+(y—a)sinip}, (26) 
where y,, y, are real phase factors which tend to zero when ak -- 0. 
When a> y > —a, the residue ¢, at « = kcosi is 


C a p| { — 4, sin 6) — (7 —78)sinip — 


= ee Gore 
sin(4, +728) 
eae i 
iz \) . ,.,/ cos(ak sin G@,)) 1? 
— (cos Fy — cos i8)(1 se oe 085%) 1c 54, sin 528< ce ane > > 
2a | 


cos(ak sin7f)} 


77 


\ 


Bt. ._ ., sin(ak sin 6,)) !? , 
< EXP lify + 54 Sin Gy sin iB —— exp ly, | » 
2| sin(ak sin7f)} 


x exptkixcosiB+(y+a)sinif}. (27) 
Finally for y < —a there is no pole at « = kcosif. Thus the wave height 
at large distances from the barrier (and from the edge of the geometrical 
shadow) is determined by the incident wave and equations (25), (26) and (27). 


5. Discussion 
We may consider first the form of the solution when the barriers are 
close together (ak->0). Then, for y > a, the expression (26) becomes 
4cosif oo = ee 
——— cos 34, cos 518 exp ik{x cos 1B + (vy —a)sin 1B}, (28) 
sin(9) +78) ' 
and equation (28) becomes 
) nee 
2cosiB cos 16, . : F ; ; as 
———— (sin }78 + sin }6,)exp ik{x cosif +(y+a)siniB}. (29) 
: Q 2 2 ~ 
sin(@, +78) 
As is to be expected, (28) is just the solution for the wave height in the 
shadow region behind a single semi-infinite barrier (cf. the particular case 
of normal incidence discussed by Crease (1956)). 

The second term of (29) shows perhaps rather more clearly than (27) 
the effect of rotation in modifying the ‘channel’ mode which exists when 
w = (0 (see below), whilst the first term gives the additional Kelvin wave 
which arises purely through rotation of the system whether there be one 
barrier or two. 
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When the rate of rotation becomes small (8—-), the residue given by 


(26) approaches zero whilst (27) becomes 
PI 


: lak iz \) 
a exp| — 4, sin #, + (1 cos Hy) y + log sak) | 
i _ ad 


(sin ak sin §,) ]!2 : ; 
< | ——_—__—___ expt, expikx, (30) 
ak sin 6, 


with 8 = Oin y%,. ‘This result corresponds to the usual channel mode without 





rotation, and has its counterpart in the acoustic and electromagnetic problem. 
Figure 2 shows for comparison the amplitude of the waves transmitted 





into the channel when cosh § 1 (w = 0) and when cosh f = 2 

(w/o = 1\3) for three different widths of the channel. ‘The latter value 
x x 1% 

Figure 2. Amplitude of ¢, on y —a for cosh f = 2 (solid lines), and cosh f = 1 


(dashed lines). ‘The amplitude behind a single semi-infinite barrier is shown for 
comparison (short dashed line). 


of cosh 8 corresponds to the , tide at the entrance to the North Sea (for 
2w we use the local value of the Coriolis parameter f). ‘The appropriate 
value of 2ak in this case (with 2a = 400 km, 4 = 109 m) is approximately 
0-297. ‘The figure also shows the amplitude of the waves behind a single 
barrier as a function of 4, when cosh 8 = 2. 

It may be seen that for the North Sea the amplitude differs little from 
that of a single barrier except for small angles of incidence. When 
6, = 42, the amplification of the channel wave is slightly greater than 2. 
This may be compared with the observed tide in the approaches to the North 
Sea and the tide at points on the Scottish coast. ‘lhe average amplitude 
of the tide at Rockall and in the Shetlands is approximately 2-1 ft., and at 
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points down the Scottish and North-East English coasts it is 4:3 ft.; so 
the observed amplification equals 2-3 approximately, in fair agreement with 
the calculated value. 


EXAMPLE II. CHANNEL WITH ONE SEMI-INFINITE AND 
ONE INFINITE BARRIER 


We now investigate a model in which long waves approach a semi- 
infinite barrier in front of an infinite barrier. ‘This and the previous model, 
together with that of a single semi-infinite barrier (Crease 1956), form 
three limiting cases for studying the additional effects that may arise at 
the edges of barriers in rotating systems. Again Heins (1956) has described 
the results for a similar problem in acoustic diffraction. ‘his corresponds 
to the special case w = 0. Much of the detail in this part is the same as 
for Example I and in general only the outline of the mathematical argument 
is given. 


6. DERIVATION OF THE INTEGRAL EQUATION 
Let the barriers be aty = a, -0 <x < w,andaty = —a,0<x< o. 
Again { may be expressed (by the use of the boundary condition (5)) as 
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Figure 3 


a contour integral round a contour I’, (figure 3) involving a Green’s function 
to be defined explicitly below; thus 


; ] oe ae OG OG 5 
C(x,v) = ——] | < C(x 9, Vo) { cos78 — +siniB — ]G,> dx, + 
; cosip Fog, ‘ ; CVy OXy Ju, =a 
{ z a C or C | 
+ | <[C]{ cos7B — +siniB — ]G,> dx, \+ 
Lie CVo OX Ju a 
yf. 2 ae . 
1G, =— =—} ds, (31) 
ONg ony 


where the arc in the last integral is the semi-circle shown in figure 3. 
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Now let G(x, y; %, Vo) be chosen so that 


re Sf . 
(cos iB— + sini )O\(u Xe 36) = 0 (32) 
OVo OX 
on y =a. With this definition it is found that the integral round the arc 
is the sum of the incident wave ¢,; and the wave ¢, reflected at y=a. If 
¢; = expik(xcos@,+ysin 4), then 
sin(4)—2 : 
= eo ad cos 6, — (y — 2a)sin 6}. (33) 
sin(4, +78) . 
‘hus a plane wave reflected at an infinite barrier in a rotating system suffers 
a change of phase. ‘The two waves together form a wave of Poincaré type 
(Proudman 1953, p. 265). 
Equation (31) now becomes 
i 1 ie 6S “9 Oo ee 0 a 
C(x, v) = — < [C]{ cos78 — + sinz8 — Gy} dx, + 


cosip 0 OVo OX eer 


u™ 





+ exp 1k(x cos 6) + y sin 9y) + 
sin(4, —7B) : 
—_——_— exp itkix cos 6, —(v—2a)sin@,'. (34 
sin(O, i iB) p ‘ 0 ( Os ( ) 
By operating with cos 78(cos 78 d/ey — siniB ¢/ex) on equation (34) aty = —a, 
v > 0, the single integral equation for [C(x , —a)] of Wiener—-Hopf type is 
found to be 


g(x) = q(x) 4 (cosis = —sinip =) : 


oy oy 

x | Floy\(cosip — +omep — Jes] dx,, (35) 

. J OVo CX V=Vo a 

where 
f(x) = 0 ree, 

a [Z (xo, 2 a)] x > 0, | 

2(x) —— Q xc > 0, 
q(x) = () + eG QO, f (36) 


= 2kcosif sin(4, —78)sin(2ak sin 4,) x 
x exp 7k(x cos 4, + 2a sin Gy) x> 0, | 
and g(x) is defined by equation (35) for x < 0. 
7. "THE GREEN’S FUNCTION G,(x, V3 Xp; Vo) 

‘The Green’s function satisfying (32) can be derived from the free space 
Green’s function of (6) by expressing it as a contour integral (cf. Watson 
1944, p. 178) of the form 

. sain i a 

LIAO R[ (x — x9)? + (v—Yo)?]!2} = -— — | x 


oO 4 4a 
Li 1x 1 


< exp thi (x —x,)cos a +(v—%)sinu] du, (37) 


where v—x) = Reos%, yv—vyp = Rsing’. The contour of integration C 


may be shifted by + 7/2 in the direction of the real axis. 
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It is then clear that G, may be represented by 


" f : p 
G,(*, 93 Xu Yo) = — | «expik[(x—x9)cos u + (y—yo)sin u] + 
FEX 
sin(u +7) 


sin(u — 7B) 


TT 
5 

exp ik[(x—x,)cosu + (v+¥9—2a)sinu]} du. 
i i } 


(38) 
The second term in the integral represents a source at the image point 
(x9, 2a—‘y) with a suitable phase shift to satisfy the boundary condition. 
The Fourier transform of G, regular for |.7{x}, < %{k} is found to be 
; sj 
n in (uo + #8) : 
sin(U») — 718) 


x exp i(k? — «?)!2!y + yo — 2a \, (39) 


K,(2,9,9) = {exp i(k®— a2)! y— yy 


a aia 
2(k? = an} 2 


where sinu, = —(k®—«?)!2/k and that branch of (k?—«?)!? is chosen 
which approaches k as « approaches zero. 


8. SOLUTION OF THE INTEGRAL EQUATION 
Taking the Fourier transform of (35) yields formally 
. . oe Seer — eee 
G(x) = O(«)+ F(a)( cos iB = ~insiniB)(cosia 2 ~ insinip) x 
Vo oy 
x FG 9, Velyin'y, » ~a 
= O(a) — F(«)L(«) (40) 
say, in notation similar to that used for Example I. 
The argument now is essentially as for Example I (there is only one 
equation to solve, however). L(«) is decomposed into factors L~(x)/L*(«) 
with appropriate regions of regularity, and we find that 


F(a) = O(a)L+(kcos )/L-(x). 





Specifically, 
Ibhcos? 
O(a) = Pt sin(9, —7B)sin(2ak sin 6, )exp(zka sin 4,), ] 
be i(a—kcos 4) | ” 
2. MBcnets 
L(x) = Ee = ‘) sin 2a(k® — x”)! exp 2ia(k? — x7)!?, 
and 








: 4 252 . . 
“Ti 2) oo) 
iia n?x ni nit (42) 
k-—«@ 


1 ; iis 1 = 
—— = (a+ kcosiB)exp| (k? — a”)! tan (=) ~ ral) | x | 


2 4a*k? 2iaxw. 2iau | 
7 ctl MO _— | 
Tt] (Fee) leSe) 
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where 


2iaz. ak 
Xo(%) ——|{ 1—y+log — }. 


7 7 


9. THE WAVE HEIGHT € 
Equation (34) may be expressed in terms of known Fourier transforms 


iS 
l 
&(x:4 ) ; . 
27 cosip 
| F(a)( cos 18 — —txsin iB) K(a390) | dz + 
P CVo Yo a 
+exprk(x cos 6) + y sin 4) 4 
sin(§, —1B ; : J 
—-—_—~ exp 1k(x cos 6, — (y — 2a)sin 9), (43) 
sin(6, +16) 
where Jik' <e <O(or < Sikcos6,'! < 0). The wave motion at a 


distance from the edge of the semi-infinite barrier is easily found by 
evaluating the residues of the integral in (43) at its poles. 

First consider the pole at x= kcos4. This will give rise to the 
‘geometric’ terms ¢, as we have seen in Example I. When these are 
combined with the incident and reflected waves, we obtain the following 
results: 


(a) x>U, a>y> a, 


ib) #2 >), a>yv> 0, 
sin(§, —78) 
C, = expik(x cos 0) + sin 6)) + ——"—~ x 
, sin(4, +78) 


x exp thk(x cos 4) — (y + 2a)sin 4) ; 


tC) £<OU, ao y> — a, 
J, sin(@, —18 
ic exp tk(x cos 4, + ysin 4,) + en eo x 
: sin(@, +78) 
x exp tk(x cos 4, — (vy — 2a)sin 4). 
If) < 2ak < z, the only other pole leading to undamped waves at a distance 


from the edge of the barrier is at x = kcos7z8. Again, when x < 0, there 
is no contribution, and the residue is also zero when x > 0, -a>y> —o. 
For x > 0, a > y > —a, the residue ¢, at this pole is 


2cos78 sin 4. sin(2ak sin 8,) ]'? /Fsin(2ak siniB) ]!? 
sin(@, +78) 2ak sin 4, 2ak siniB 


a : 
2iak “0 E 17 
x EXP - | (cos — cos a1 —y+log zat) - 


#) sin 6) — (7—78)sin 7B 4 ts Jexp tk(xcosiP +(y+a)siniB), (44) 


where ¢, is the same as in Example I, except that 2a replaces a, 
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10. Discussion or EXAMPLE I] 

There are certain interesting differences between the wave amplitudes 
in the channel found in this and the previous example. 

First, as 4, approaches zero in (4+), the amplitude approaches zero. 
Mathematically, the explanation lies in the expression for the reflected 
wave given by (33). This wave suffers a phase change relative to the 
incident wave, which approaches 180° as 6) becomes small. The effects 
of the incident and reflected waves thus cancel one another. However, 
the asymptotic form of (44) as 4) approaches zero will be a good approxi- 
mation only at increasing distances down the barrier. Physically speaking, 
the particles in a plane wave in a rotating system have transverse accelerations 
which are balanced by the Coriolis force due to motion in the direction of 
propagation. When 4, approaches zero, these transverse accelerations 
would be in a plane almost perpendicular to the barrier. As there can be 
no motion through the barrier, these accelerations must be annulled by 
those from a reflected wave almost 180° out of phase with the incident 
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Figure 4. Amplitude of ¢, on 4 —a for cosh f = 2 (solid lines) and cosh f = 1 


(dashed line). 


Second, the term which in Example I was found to arise essentially 
from the rotational effects at a single semi-infinite barrier (cf. equation (29)) 
is seen to be absent from equation (44). ‘This again may be understood by 
considering the rotational effects of two plane waves approaching a single 
semi-infinite barrier, one at an angle 6) and the other at — 6) with a phase 
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shift of sin(#,—78)/sin(@)+78). It is found that the waves propagated 
down the barrier as a result of these two incident waves annul one another. 

Figure 4 gives the amplitude at the barrier y = —a, x > 0 for cosh B = 2 
corresponding approximately to the M, tide in the latitude of the North Sea. 
‘The value of 2ak for this sea is 0-297. It is seen that the amplitude of a 
normally incident wave is about 3-0, which is in excess of the observed value 
of 2-3 quoted for Example I. 
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The thermal boundary layer on a non-isothermal 
surface with non-uniform free stream velocity 
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N.A.C.A. Lewis Flight Propulsion Laboratory, Cleveland, Ohio 


(Received 31 December 1957) 


SUMMARY 

A formally exact solution for the thermal boundary layer on 
a non-isothermal surface subjected to non-uniform free stream 
velocity is presented in the form of a series. It is demonstrated 
that the solution can be recast in terms of universal functions, 
which are independent of the wall temperature data of particular 
problems, and which depend only on a single parameter character- 
izing the variation of the free stream velocity. 


INTRODUCTION 

A new method for computing steady laminar incompressible boundary 
layer flows under conditions of non-uniform free stream velocity has recently 
been formulated by Gortler (1957). Provided that the free stream velocity 
can be written in terms of rather general series (or polynomials), he is able 
to give a formally exact solution for the boundary layer velocities. Further, 
he shows that the solution can be written in terms of universal functions, 
which depend only on a single parameter of the free stream flow. 

Our interest here is in the thermal boundary layer, and we consider 
the general situation of a non-isothermal surface with non-uniform free 
stream velocity. Gdortler’s work will constitute a point of departure for 
the present study, and we will be able to make direct use of his results. 

For situations where variations in both the surface temperature and 
free stream velocity are expressed as certain rather general series (with 
arbitrary coefficients), we are able to give a formally exact solution for the 
boundary layer temperature distribution. Our final results will involve 
a new set of universal functions. When these have been determined, the 
computation of the heat transfer and of other thermal quantities is a simple 
algebraic process. 

Exact solutions of the boundary layer energy equation for non-isothermal 
surfaces have previously been found only for free stream velocity variations 
in which U « x”. ‘These are the so-called similar solutions. ‘The present 
formulation permits the consideration of much more general variations 
of free stream velocity. 

The development is carried out for steady laminar flow with constant 
fluid properties. Later, modifications for variable fluid properties will be 
noted. 
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‘THE ENERGY EQUATION AND ITS TRANSFORMATION 
The differential equation expressing conservation of energy for steady 
laminar non-dissipative flow in a boundary layer is, in the usual notation 
(x being the thermal diffusivity), 


.=— fVU— = a—. (1) 


‘The boundary conditions for the temperature are assigned as 
T(x, 0) = T,,(x), (2a) 
T(x,y)>T,, as yv/(Re)> a, (2b) 


where 7, is a constant. ‘The terms describing the frictional heating and 
compression work have been omitted from the energy equation, thus 
restricting the problem to free stream velocities for which the adiabatic 
temperature rise is much smaller than the temperature difference between 
surface and stream. 

An approach to the solution of equation (1) cannot be made without a 
knowledge of the velocity components u and v. A consequence of the 
constancy of the fluid properties is that a complete solution for the velocity 
may be made without recourse to the temperature. This same desirable 
independence of the velocity is also found for a special group of variations 
in fluid properties to be noted later. So, in relation to equation (1), the 
velocity components u and wv can be regarded as known a priori. Since our 
goal is to attack as broad a problem as possible, we will select the most 
general velocity solutions presently available: namely, those of Gortler 
(1957). 

Following Gortler, we introduce new independent variables € and 7 
and a new dependent variable F by the relations 


- ( a ye 
é = : | U dx, n = yUs ay | Udx> , (3) 


c 


F = p/vy (2£), (4) 
where u is the stream function. In terms of these variables, the velocities u 
and wv are 
u= UF, v= —U(2€)-'*[F+ 2éF, + (B—1)yF_], (5) 
where the subscripts denote differentiation and 8, termed the principal 
function by Gortler, is given by 


p= | Ude. (6) 


For the temperature problem, we introduce a dimensionless variable © 
by the definition 
r-7, . ~Te (7) 





 —. a 
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Using the new independent and dependent variables, the energy 
equation (1) and the boundary conditions (2) may be written as 


T° 
Pr)0,, + (F+2€F.)0, = 2éF | Of —*.) +0, |, 8 
(Pr)10,,, + (F+2€F,)®, = 26F,] 0( 7.) +0, (8) 
®(E, 0) = |, O(E, n) > as n> ®, (9) 
where 7’ is an abbreviation for d7,,,/dx and Pr is the Prandtl number. 


All the given data concerning wall temperature and free stream velocity 

in any particular application appear explicitly only as the single coefficient 

vT’ UAT, which we will call the principal thermal function and denote 

by A. It is interesting to note that A arises in the temperature problem in 
> 


a manner similar to that in which £ arises in the velocity problem. Both A 
and § will play important roles in the solution of the temperature problem. 


GORTLER’S SOLUTIC’'S FOR THE VELOCITY 
Asa final prelude to the solution of equation (8), the required information 
for the velocity variable F must be given. We choose to study situations 
where the free stream velocity may be represented by the convergent 
series 
U(x) = x {59 +5, 2+! + sq e8Mt +...) —l<m< ao, 5 #9. (10) 
Included here are two very interesting physical situations. First, the case 
m =() represents flows with a forward cuspidal point at v = 0; and in 
particular, the flow over a flat plate is given by s, = s, = ... = 0. Secondly, 
the case m = 1 represents flows about a symmetrical body* with a forward 
stagnation point. 
The form of the principal function 8, which corresponds to the free 
stream in equation (10) is 
B = By t+ B, E+ Bo? +..., (11) 
where 8, = 2m/(m+1) and the other £,, depend on the coefficients s,, of 
equation (10). Further, the relationship connecting the coordinate a 
with & is found from equations (3) and (10) to be 
tm hl > Sn yim +), (12) 
v(m+1) Ay n+1 


Finally, the solution for F is given by Gortler as 
F(é,n) = X F,(n)é", (13) 


where in turn the F,, are written as the following linear combinations of 
universal functions 


F, —* By fis ) 

F, = Bifir+ Bohs» | 

F; thin ' 3 B, fi2 + Bsfs, (14) 
Fy = Bi fia + PF Poh + Pr Bs his 33 foe aif 


* As pointed out by Gértler, the restriction of symmetry can be removed by con- 
sideration of an alternate series for the free stream velocity. 
X2 
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‘The function F, is determined from the Falkner—Skan differential equation 
containing By as a parameter, while f,, f\;, fo, etce., are found from the 
solution of linear inhomogeneous ordinary differential equations with 
homogeneous boundary conditions. ‘The important fact is that the 
equations for the functions f do not contain f,, By,.... Only By enters 
indirectly through the appearance of Fy and its derivatives. So, for a 
selected fy, all the functions f as well as Fy can be computed once and for 
all, without further reference to the additional special data (f,, Bs, ...) of 
the particular problem. It is in this sense that the functions f are universal. 


SERIES SOLUTION OF THE ENERGY EQUATION 
We seek a solution of the energy equation (8) in the form of a series as 
follows: 
® = p2 4 (n)é". (15) 


We will assume that the principal thermal function A may also be expressed 
as a series: 


A = (16) 


Ms 
‘as 


0 
It will be shown in a later section that the dependence of the wall temperature 
on x which corresponds to this form for A is a rather general series expansion. 
Further, we will show how to determine the coefficients 4,, for prescribed 
variations of wall temperature and of free stream velocity. 

Introducing the series (13), (15) and (16) into the energy equation 
and grouping terms having common powers of €, we find a system of ordinary 
differential equations for the 6, : 

(Pr) 100+ F,0=0, (0) =1, (x2) =9, (17) 
(Pr)*0, + F,0,—2nF,0@,,=¢,, 6,(0)=0,(0)=0, (= I, 2,...), (18) 


where 
. . vv . / " i e. wv 
$. = 2 (2a—p)F,¢,— (+ 1)F, 4G, .,1+2 2 An» ( » F;9, i) 
pl p=0 j—0 
If we introduce the abbreviation 
M,,(y) = (Pr)-1y" + Fay’ —2nF;, y (19) 
then the first few of equations (18) are 
M,(0,) = —3F, 6, + 2A» Fo 9%, (20a) 
M,(0,) = 2F, 0, —3F, 0, —5F 6, + 2A9(F 0) + Fy 0) + 2A, Fy 9%, (20 b) 


M5(45) 4F, My =3 F, 0, ! oF, 0; 7 SF 0, * 7F, 6, U 
1 2A, (F,Ay + FY A, + F18,) +24 (F% A, + F8,) + 2ry F 8p, 


(20) 
Equations (17) and (18) constitute a recursive system for computation of 
the 6,. All of equations (18) are linear. The boundary condition at the 


wall (7 = 0) is satisfied by 9 alone, 
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It is interesting to inspect the function ¢,, which contains the input 
data for the computation of #,. It is seen that the particular data embodied 
in By, -) 8B,» Apgy +) A,,-3, and the Prandtl number, must be specified 
before equation (18) can be solved for 6,. It is desirable to rephrase our 
problem in terms of functions which are independent of the explicit data 
of a particular problem. Such a course is followed below, where we will 
find universal functions which are completely independent of particular 
wall temperature conditions and which depend on the free stream velocity 
data only inasmuch as they depend on {,. 


The universal functions 
Our aim is to reduce the functions @,, to a linear combination of universal 
functions. First, we consider 6, and its differential equation (20a). The 
velocity functions appearing on the right-hand side may be replaced by 
universal velocity functions (14) to give 
M,(0,) = —B,(3f, 9) +Ao(2F% 9). 
Inspection of this equation suggests that 6, be written as 
9, = B, By +Ag Ly. (21) 
where B, and Ly are functions of 7 which satisfy the differential equations 
M,(B;) = 3f; 9, B,(0) = B,( 0) = 0, 


(22) 

Mi(L,) = 2F,4,, L,(0) = L,( 20) = 0. 
Since Fy, f,, fi, etc., depend only on 8p, it follows from (22) that the same 
is true of B, and Ly. These functions do not depend upon particular 


surface temperature data. ‘Thus, they are universal in the same sense as 
are the velocity functions. 

Next, consider 9, and its differential equation (20b). When the right 
side is evaluated using the universal velocity functions (14) and equation (21) 
for 6,, we find that the factors B?, B,, 8, A,, A=, and A, appear. So we write 

6, = B? By, + By Bo + By Ag Z1,9 +A2 Log +A, Li. (23) 
The differential equations for B,,, Bs, etc., are found from the usual method 
of comparing coefficients to be 
M,(B,,) = 2f, By —3f, By — 5h 9), 
M,(B.) = —5f,9,, 
M(Z,0) = 26, (99+ Lo) —3f, L,+2F, B,, } (24) 
M,(Loo) a 2F, Lo 
M<(L,) = 2F,0,, 
where 
B,,(0) = Byy( 0) = Z, (0) = Z,o( o)=..= 0. 
Again, we have universal functions in the sense noted above. 

In this fashion we can continue to recast each 9, as a linear combination 

of universal functions. Further listing of the universal functions and their 
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corresponding differential equations is given in Appendix A. It is not at 
all surprising that there are considerably more universal functions involved 
in our present velocity-temperature problem than there were for the velocity 
problem alone (compare equations (21), (23) and Appendix A with equation 
(14)). 

It is interesting to look at some of the properties of the universal functions. 
An inspection shows that the differential equation for any one of the 
functions B does not include any of the functions L and Z as input data. 
The functions B may therefore be computed independently of the 
functions L and Z (the converse is not true). Using the results of Sparrow 
(1958), it may be observed that the functions B are simply the universal 


functions for an isothermal wall subjected to the free stream velocity of 


equation (10). ‘The functions Z and Z are associated with the deviations 
of the wall temperature from a uniform level. ‘This splitting of the problem 
is, of Course, associated with the linearity of the energy equation. 

‘The symbolism used for the universal functions follows a rational 
path first set down by Gortler. It permits an easy way of identifying any 
universal function with its coefficient as they appear in the linear aa’ 
which compose the 6. For example, B, is multiplied by f3, Lo, by Apo 
and A,, B,, twice by f). 


Series representation of wall temperature variation 

We now investigate the wall temperature variation corresponding to the 
assumed series representation for the principal thermal function A. We 
had written 


4 hag 


A = = oS A é". 1¢ 
a es at i. ne ( >) 


For the velocity variation a the can idl between € and w Is given 
by (12). Using these facts, — (16) may be rewritten as 





w+1 a . Pp \ 
— = = yl ee yn 1) : I + Sn gine 1) . (25) 
T,—T, v(m+1) ,“on+1 
It is easily seen that si series for oT, which satisfies equation (25) has the 
form 








T,.-T, = Sa,x — (a, #0). (26) 


t 


For the two important special cases of flow with a forward cuspidal point 
(m = 0) and of flow around a symmetrical body with stagnation point (m = 1) 
equation (25) reduces to 


T,-T., = > 4," (a,#0, m=O), (26a) 
0 

T,-T, = > a,x" (440, m=1). (26b) 
1=0 


For particular problems, the wall temperature and free stream velocity 
ands, will be specified. It is necessary 


variations will usually be given, 1.e. a 


” 
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to determine the A,, from these. ‘The needed relationship may be evaluated 
from (25) and (26). The first two of these are written below; a more 
complete listing is given in Appendix B: 


v(m+1 = 
ro = = da,, (27 a) 
So 4 
2 1)2 ; 2 
-_ a 2a, <3 |. (27 b) 
Sj) Ao So 4 


It is clear that the determination of the A, from the given data is simply 
an arithmetic process. 


Heat transfer 
When the wall temperature is specified, the quantity of greatest practical 
interest is the heat transfer. From Fourier’s law, the local heat flux g at 


the wall is 
q=- (: =) ‘ (28) 


( ry 


In terms of the variables of this analysis, the expression for g becomes 


qx 2 | v A’ n a 
Ss & (O)E 29 

(T,- T’,.)k Ra ( ) 
where Re is the Reynolds number Ux,v, and 6 (0) is an abbreviation for 
[dé,,/dy]y. ‘The 0° (0) may then be replaced by the universal functions in 


21), (23) and Appendix A. 
) ) Pp 


Variable fluid properties 

Both Gortler’s analysis and the one given here apply, with no essential 
change, to fluids with variable properties under the following circumstances : 
(1) pu = constant, (2) pk = constant (k = thermal conductivity), 
(3) ¢, = constant. ‘Then, using Howarth’s transformation 


Y = | (p/p...) ay, (30) 


and replacing y by Y in the equation (3) defining », the entire formulation 
of the solution remains as before. 


CONCLUDING REMARKS 

We have given a formally exact solution for the temperature problem 
under rather general variations of wall temperature and of free stream velocity. 
As a consequence of the introduction of universal functions, application of 
the results to particular situations becomes an arithmetical procedure. ‘The 
practical utility of the development rests upon the availability of the universal 
functions. At present, only the functions B for Pr = 1 which correspond 
to the important cases of m = 0 and m = 1 are available (Sparrow 1958). 
Efforts are being made to find high speed, large memory computing equipment 
on which to evaluate a larger group of universal functions. 
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ApPENDIX A, UNIVERSAL FUNCTIONS AND THEIR DIFFERENTIAL EQUATIONS 
‘The universal functions are introduced by the following equations: 
6, = B, By +Ap Lo, 
0, = 8? B,, +B. B.+B, MZ 9 tr2 Log +A, Ly, 
6, = B3 Byxy + B, Bp Bys + By By + B2Xy Zrn.9 + Bi A2 Z,50 + Bry Za 4 
+ Ba Ay Zo9 + Ag Le +Ag Ay Lo, +A3 Logo 
1 Bz Bry + By Bs Big + Bg Boo + By By + BEAy Zino + 
AS 411,00 ' BrA TAT + By 32Ay Z12 2,0 T By AZ, 2 : 
t By Ag Ay 21.01 + Bi AG 21.000 + Be A3 20.00 + BeAr 201 4 
+ Bs Ag Z3.9 + Ag Lg + Ag Ag Log +A2 Ay Loo + Az Ly, + At Loooo- 


} B34 i 
y BY Buy 


oD 


Using the W_,(y) operator defined by (19), the differential equations for the 


universal functions are as follows. 
Casen = 1. See (22). 
Casen=2. See (24). 


Case n = 3. 
M;(Byy1) = 4f; By — 3f; By, y 2B, fi By 7hr11 9% 
M,(By,) = 4f; By — 3f, By + 2f, By — 5fz By — 7h % 
M,(B3) = —7f,9, 
M3(2119) = 2F, By, + 2f; (By, + 2210) — 3f, 21,0 + 2f;,(% + Lo) — afr Ly 
) = 2Fo 219+ Of, (Ly + 2h) —3f, L 
) = 2F, By + 2f; (9% + 2L,) —3f, Li, 
) = 2F) Be + 2f; (09+ Ly) — 5fo Ly, 
) = 2F'0,, Me(Lo1) = 2F(Lot+L;), Ms(Looo) = 2F/ Loo: 


t 
00” 


Case n = 4. 
M,(Bii) . 6f; Bi 3f; Bes y 4h; By, Shi B, oe 1 B, Thi B 
a Of1111 9% 
M,(By2) = Of; By.— 3h; By, t 4h By— Sf By + 2fy2 B,- 7h By i 
oe 4f. By, — Sf By, ~ Hrr2 Oo» 
13) = of; B, 3f; By fis 0, y 2f; B, 2 7h B;, 
>) = 4f5 Bz —5f2 Bo — fos 9, 
)= 2 = Y, Ais 
,= = 2F, Bu T 2f (Bir + 32310) —5f; 


1,(B 
mt 
1,( 


ee “111.0 ob et 2210) 


a Sf 2 Z,0* Aino + Lo i Thin Lo 
M (211.00) : TF, Bis t 2f,(Z 41.07 321.00) = 3f, Z. “1,00 7 stg T ZL) 7 
—Sful 
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M4111) - 2Fy By, 4 2, (B, ' 34,1) ‘i Sf, Zi. ; 2f11(% i 2L,)- Si Ly, 
My(Zy20) oe 2F, By» Ls 2f; (Be r 322) - 3h; Z39 r 2fio(4 + Ly) r Tha Ly 7 
+ 2f,(B, ss 22,9) i Shs Zi. 
M,(Z, 2) = 2F) By + 2f; (9 + 3L2) — 3f, Ly, 
M,(2Z,.1) = 2F,(Z10+ 21,1) + fq (Lo + Ly, + 3291) — 3h Bows 
M,(Z, 000) = 2F 21.00 + 2h (Loo + 3Lce0) — 3h; Lewes 
M, (Zo) = 2F, Z20 + 2fe(Lo + 2Ly9) — She Lows 
M,(Ze,) = 2Fi Be + 2f,(0,+2L,)— 5fgL;, 
M,(Zoo) = 2F, Bs +2f,(0+Lo)— fale,  M,(L3) = 2Fi 0%, 
M,(Lo2) = 2F(Le+Ly),  My(Loo1) = 2Fi(Loo + Ler), 
M, (Ly) = 2F,L1  My(Looo) = 2F Looe: 


APPENDIX B. DETERMINATION OF ),, FOR GIVEN S, AND 4, 


¢ n+l 
sien hy =a al a7 | : 
v(m+1) 
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On the breaking of water waves of finite amplitude 
on a sloping beach 


By H. P. GREENSPAN 
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SUMMARY 

In a recent paper Carrier & Greenspan (1958) showed that, 
within the framework of the non-linear shallow-water theory, there 
exist waves which do not break as they climb a sloping beach. ‘lhe 
formation of a shock or bore is dependent on a variety of factors 
(wave shape, particle velocity, etc.) and, as yet, no general criteria 
for breaking have been found. In this paper, we consider waves 
which propagate shoreward into quiescent water; it is shown that 
any compressive wave (a wave of positive amplitude) which has 
a non-zero slope at the wave-front eventually breaks before reaching 
the coastline. In fact, an explicit relation is obtained between the 
initial conditions and the position where breaking occurs. 


‘The conservation equations of mass and momentum of the non-linear 
shallow-water theory are 


[v*(* : al | ag = — ye, (1) 
and Use +u*ut, = —Z" Ne (2) 


where the symbols »*, A*, and x* are defined in figure 1, v* is velocity, 
t* is time, and g* is the gravitational acceleration. ‘The asterisks denote 
dimensional quantities. Let the depth be given by 


h* = a(l* —x*), 
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Figure 1. Fluid with a fixed boundary and a free surface. 
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and v* = (xg*/*)'*. With the substitution of these dimensionless variables 


and let x = x*//*, 1 = y*/al*, v = v*/v*, t = t*/t*, where t* = (/*/ag*)!? 
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the basic equations become 
[e(l—x+y)], = —%% (3) 
and U+UU, = —,. (4) 
‘hese hyperbolic equations can be rewritten in a form in which the 
characteristics «, 8 play the role of independent variables and v, x, t and 
the wave velocity c = (1—x+ 7)!? are unknown functions of «, 8. The 


characteristic equations which arise are 


X_ = (vteyt,, (5) 

(v+2c+t), = 0, (6) 

y= (v—e)t,, (7) 

and (v—2c+t), =0. (8) 


‘he reader is referred to Stoker (1948) for a complete derivation of these 
equations. ‘I'wo families of characteristics are described. ‘The quantity 
v+2c+t is constant along any characteristic « = const. which propagates 
shoreward; the slope of this characteristic at any point (x,/) on it is 
dx/dt =v+c. The quantity v—2c+t? is constant along any characteristic 
8 = const. propagating seaward; the slope of the characteristic at any 
point on it is dx/dt = v—e. 

Consider, then, waves which are propagating shoreward into quiescent 
water, and which at time ¢=0 are given by 7 = 0 for 0< x < 1, and 
7 = f(x) for x < 0, where f(x) is a known function and f(0)=0. The 
characteristic forming the wave-front can be determined explicitly. Since 
the wave is propagating into a region of rest, the particle velocity wv is 
identically zero at any point on this characteristic. ‘Therefore the slope 
of this curve is given by dx/dt=c. ‘The quantity v+2c+t is constant 
along this characteristic so that v+2c+t = 2c+t = 2, since c= 1 at x=0 
and t= 0. ‘This implies that the wave-front characteristic is the curve given 
by dx/dt = 1— 3$t, with x = 0 at t = 0; that is, the branch of the parabola 
x = t—4t? for which ¢ < 2. 

The wave velocity at the wave-front is 

e = (1—x)¥? = 1-32. (9) 

At zero time the wave-front is located at the origin of the coordinate system 

fixed in the fluid. At a subsequent time ¢ the wave-front has moved a 

distance c dt = t—}t? from the fixed coordinate system. If € measures 
4/0 


distances from the moving wave-front (see figure 2), then 


=E+ | cdt =&+t—}#*. (10) 
0 


& 


‘The origin of the new coordinate € is the position of the wave-front. If we 
now introduce (10) into the basic equations (3) and (4), we find that 


[w(1—€-t+40?+)], = —m+(1—3t)m, 
= u(1—£—t+}i?+y)+u(—1+7,), (11) 
and u,—(1— 3t)u,+uu, = —7,. (12) 


In the moving coordinate system, 7 = 0, 4, = 0, u, = 0 and uy, = 0 at € = 0. 
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If we set € = 0 in equations (11) and (12), we obtain the result 
He = (1— 3t)u, at € = 0. (13) 
By differentiating (11) and (12) with respect to ¢, we find that 
la = Mg(1—3t)— Ju, at = 0. (14) 
Similarly, by differentiating with respect to €, we find that 
(1 — $¢)—2u,+3(1—$)eg = —n¢ af = 0. (15) 


By eliminating »,, from (15) by means of (14), there results an equation 
for u,(0, t) alone: 


(1—}t)ug = gu,—3(1— 3t)ue at € = 0. (16) 
‘The corresponding equation for 7,(0, t) is 
(1— 3t)ng = dng 20%. (17) 








ail Nn 
A A 
y o> es 
A te 
+ —__*+—_——__ > 
< - — > 
rt 2 
J cdt=t-t 
J; Fe 





Figure 2. Location of the coordinate system fixed with respect to the wave-front. 


Thus, if 7,(0,0) < 0, the wave is compressive at the wave-front, and 
(9,0) < 0 or the wave-front steepens. If the wave is initially compressive, 
with a non-zero slope at the wave-front, the wave-front steepens as it advances 
on the shoreline. On the other hand if 7,(0,0) = 0, then 7,(0,t) = 0 
for all subsequent times. Such waves cannot begin to break or form a 
bore at the wave-front. Indeed, they may or may not break at some interior 
point. If we consider rarefaction waves for which 0 < 7,(0,0) < 3, we 
find that such waves actually steepen at the wave-front in contrast to rare- 
faction waves on a constant depth ocean, which always flatten out. A wave 
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for which 7,(0,0) = 5 continues to propagate shoreward with this slope 
at the wave-front. Waves for which 7,(0,0) > } flatten out as they advance, 


and the slope at the wave-front approaches }. 
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Figure 3. Breaking time as a function of initial slope. 














Figure 4. Breaking position as a function of initial slope. 


Equations (16) and (17) are non-linear first-order differential equations 
for the functions u,(0,t) and 7,(0,¢) with the boundary conditions that 
u.(0,0) and »,(0,0) are specified. The solutions of these differential 
equations are 

u(0,t) = 1/{(2—t)[1 — A221 — 3t)8?]}, (18) 
and n (0, t) = 1/{2[1 — A??(1 — $2) }}, (19) 





334 H., P. Greenspan 


where 


u.(0,0)—172  [»,(0,0)—472 
u-(0, 0) ~L 7,(0,0) 
If 7.(0, 0) m where m > 0), then A is larger than one. A wave satisfying 


this condition breaks when the slope at the wave-front becomes infinite. 
From (19) it is seen that this occurs when 


2m “78 
ee eee oe ee Oe a ee 
| F ; a J 


or equivaiently at x = t—}f < 1. 

‘Therefore waves which are compressive in the region adjacent to the 
wave-front and propagate shoreward with a discontinuity in slope eventually 
break before reaching the coastline. Values of t and x for which breaking 
occurs are plotted against m in figures 3 and 4. 

Although general criteria for breaking are still lacking, the breaking 
point of a compressive wave with small radius of curvature near the 
wave-front can be accurately predicted. 


This work was sponsored by the Office of Naval Research under 
Contract Nonr 1866(20). 
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REVIEWS 


Fluid Models in Geophysics, edited by Rosert R. Lone. 
Washington: U.S. Government Printing Office, 1956. 162 pp. 
$2.25. 

The book consists of papers presented at the first symposium 
use of models in geophysical fluid dynamics ’’, at Johns Hopkins University 
in 1953. All deal with hydrodynamical problems with some relation 
(occasionally rather distant) to geophysics. ‘The first, by Stanley Corrsin, 
is a general introduction to dimensional analysis, with special reference to 
the non-dimensional parameters that appear in the statements and solutions 
of the problems, and must be the same in the model as in the original if the 
solutions for one are to be adapted to the other. ‘The principle is used 
regularly in aeronautics, but geophysics is still cluttered with experiments 
that are irrelevant because similarity conditions are not satisfied. ‘They are, 
I think, satisfied in all the problems considered in this book except possibly 
some that are explicitly stated to be of an exploratory character. ‘This 
chapter would have been clearer if some passages that suggest that the 
theory of dimensions depends on units had been omitted; the formal 
analogy obscures the essential difference between the applications. 

‘The remaining papers are applications of hydrodynamical methods to 
problems of meteorology and oceanography and of motions in the central 
core. The last subject has probably some connection with terrestrial 
magnetism and is considered by R. Hide. The Elsasser—Bullard theory 
of the magnetic field supposes the field maintained by a dynamo action 
in the fluid core, the energy being supplied by thermal convection. A full 
analysis has not yet been carried out, but a model based on a rotating cylinder 
of fluid with an internal source of heat gives motions that are interesting 
for their own sake, and suggests possible actual motions in the core that 
may be relevant to the secular variation. ‘I'wo types of motion are specially 
noted and an explanation is suggested by E. N. Lorenz. 

The other papers are very difficult to summarize, because they are 
already highly condensed summaries. ‘The bibliographies given should 
be most valuable for further reading. ‘They largely concern motions 
produced thermally either by long-period disturbances or by instability, 
One specially interesting suggestion, due to Rossby and mentioned by 
von Arx, is that some problems of ocean circulation can be reduced to 
laboratory conditions by a change of the form of the boundary. If this can 

be extended we might get an approach to a theory of the tides in the actual 


‘ 


“on the 


ocean. 

One point that needs amendment concerns Proudman’s theorem that 
in a rotating fluid every small steady motion must be two-dimensional. 
What is proved is dw/dz = 0, with suitable axes; it does not follow that 
w = (0) everywhere unless there is a surface across the fluid such that w = 0 
at all points of it. Usually, of course, there is such a surface, namely the 
bottom of the ocean or of the tank, and the conclusion follows. But it 
would be possible to lift the tank up, and the fluid must then move with it, 





336 Reviews 


| came upon the point in discussing certain motions of the earth that involve 
changes of direction of the axis of figure; the ellipticity of the core boundary 
produces motions in the core parallel to the axis of rotation, which can be 
comparable with the transverse motions. 

As an attempt to compensate for a very inadequate summary, perhaps 
I might be allowed to mention a few points that might be considered at 
future symposia. 

How can we reconcile the phenomenon of ground mirage with the 
present theory of thermal instability? ‘There is a variation of temperature 
of severa! degrees in the lowest centimetre of the air. If there were fixed 
surfaces above and below, or a fixed surface below and a conducting free 
surface above, this state would be stable. But with air on top it seems 
that the steady state can be maintained only if the heat is carried away by 
turbulence, and a full understanding would apparently involve a turbulent 
atmosphere with the mirage layer as a sort of boundary layer. 

In an early paper on travelling atmospheric disturbances (Phil. Mag. 7, 
1919) I took a cyclone of the type then considered by Sir Napier Shaw, 
superposed on a general circulation, and studied the rate of change of 
pressure, taking higher powers of the velocities into account. The 


surprising fact was that with ordinary values of the velocities the rate of 


motion of the cyclone would be far too fast, and actual rates seemed to 
require that the actual isobars are much more nearly circular than in the 
model adopted. I can think of several respects in which the early model 
departs from what we now know about cyclones; probably modern 
meteorologists will think of others. Still, I should like to know just what 
properties of actual cyclones permit them to travel so slowly. 

It is generally believed that the apparent secular acceleration of the Moon 
is due to tidal friction in shallow seas, though the data have never been 
as consistent as we should like. Recent work, especially by Urey, Holmberg, 
Egyed and C. A. Murray, has claimed alterations in different ways. 
A possibility that needs serious consideration is that the dissipation of 
energy in the tides may be from 3 to 7 times what has been thought. The 
known shallow seas seem unable to explain this; what can be done with 
tidal currents along the open coasts of continents? It was always difficult 
to see how so much energy got into the shallow seas at all. In my estimate 
about half the dissipation was in the Bering Sea, for which the data were 
rather rough; but now that Alaska is a region of scientific activity better 
information should be available and a better estimate should be possible. 

Finally, I venture to suggest that something should be done about the 
language next time. Most of the papers are extremely difficult to read, 
with involved constructions that recall old-fashioned German. ‘There is no 
need for this. I should certainly not claim that English physical writers are 
perfect, and some American scientists certainly could write intelligibly ; for 
instance, H. N. Russell and R. A. Daly. Why, for instance, should the word 
‘celerity’ be dug up for what we have for a century been calling wave- 
velocity? 

HAROLD JEFFREYS 






































Sound Pulses 
F. G. FRIEDLANDER 


A description, based on the theory of linear partial differential equations of 
the hyperbolic type, of the theory of sound pulses and its recent develop- 
ments. Individual chapters deal with the equations of motion, wave-fronts 
and their characteristics, geometrical acoustics and their application to 
reflection problems, and with the diffraction of pulses. 40s. net 


Nuclear Scattering 


K. B. MATHER & P. SWAN 


Much information on nuclear forces has been derived from a study of 
nuclear scattering. ‘This book gives a clear account of the physics of nuclear 
scattering, including experimental techniques, significant experimental 
results and interpretation in terms of nuclear forces and nuclear structure. 
Of interest to all nuclear physicists, practical and theoretical. 
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